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NOMENCLATURE 

2 

A  =  Heat  transfer  area,  ft 

C  =  Solute  concentration,  1 bs/( 1 b  solution) 

C  =  Specific  heat,  Btu/lb 

D  =  Diameter,  ft 

G  =  Mass  velocity,  lbs/(ft  )(hr) 

Gr  =  Grashof  number,  dimensionless 

H  =  Hold-up,  lbs 

I  =  Index  set 

K  =  Index  set 

L  =  Length,  ft 

N  =  Number  of  tubes 

2 

P  =  Pressure,  lbs/(in  ) 

P(t)  =  Covariance  of  estimate,  vector 

Pr  =  Prandtl  number,  dimensionless 

Q  =  Heat  transfer  rate,  Btu/(hr) 

Q(t)  =  Covariance  of  process  noise,  vector 

R(t)  =  Covariance  of  measurement  noise,  vector 

Re  =  Reynold's  number,  dimensionless 

T  =  Temperature,  °F 

AT  =  Temperature  difference,  °F 

U  =  Overall  heat  transfer  coefficient,  Btu/(hr)(ft  )(°F) 

U(t)  =  Control  vector 

V(t)  =  Process  noise,  vector 

3 

V  =  Vapor  volume,  ft 

V  =  Vapor  flow  rate,  lbs/mi n 


W  ■  Liquid  flow  rate,  lbs/mi n 

X  =  State  variables,  vector 

X  =  Estimate  of  state  variables,  vector 

X  =  Lockhart-Martinelli  factor 

Y  =  Calculated  observations,  vector 

Y  =  Actual  observations,  vector 
f  =  Function  of 

g  =  Acceleration  due  to  gravity,  ft/(hr  ) 

h  =  Liquid  enthalpy,  Btu/lb 

hv  =  Vapor  enthalpy,  Btu/lb 

h  =  Film  coefficient,  (Btu)(ft)/(hr) (ft  )(°F) 

k  =  Thermal  conductivity,  (Btu)(ft)/(hr)(ft  )(°F) 

p       =  Variance  of  estimate 

t       =  Time,  minutes 

u       =  Control  variable 

v       =  Process  noise 

w       =  Measurement  noise 

x       =  State  variable 

x       =  Estimate  of  state  variable 

Subscripts 

0  =  Outside  world 

1  =  First  effect 

2  =  Second  effect 

i  j      =  From  unit  j  to  unit  i 

a       =  Before 

c       =  Condensate 

xi 


s  =  Sensible  heating  zone 

t  =  Tube 

B  =  Boiling  zone 

w  =  Wall  conditions 

F  =  Feed 

i/j  =  At  "i"  given  conditions  at  "j" 

f  =  Film  conditions  or  final  condition 

Superscripts 

i  =  Inside 

o  =  Outside 

v  =  Vapor 

*  =  Optimal ity 

Greek  Letters 

a  =  Point  constraint  multiplier 

3  =  Inequality  constraint  multiplier 

0  =  Estimated  parameters,  vector 

6  =  Estimated  parameter 

p  =  Density,  lbs/ (ft3) 

A  =  Lagrange  multiplier  or  latent  heat  of  vaporization, 

Btu/lb 

u  =  Viscosity,  lbs/(ft)(hr) 
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The  application  of  optimal  control  theory  to  a  chemical  engi- 
neering problem  is  investigated  by  the  development  and  implementation 
of  a  control  policy  for  the  minimum  time  start-up  of  a  double  effect 

evaporator. 

The  particular  evaporator  on  which  experimental  runs  were  made 
was  a  laboratory  scale  double  effect  evaporator  with  reverse  feed.  It 
was  completely  instrumented  for  control  by  the  installation  of  orifices 
for  measuring  flow  rates,  thermocouples  for  measuring  temperatures, 
pressure  taps  for  measuring  pressures  and  hold-ups,  and  pneumatic  con- 
trol valves  for  manipulating  flow  rates.  Transducing  and  controlling 
instruments  were  installed.  In  order  to  do  on-line  computerized  data 
logging  and  control,  interfacing  of  the  process  with  the  IBM  370/165 
computer  on  campus  was  provided  by  an  existing  IBM  1070  interface.  A 
part  of  the  existing  user  circuitry  associated  with  the  interface  had 
to  be  rewired  and  modified  to  function  appropriately. 

The  dynamic  model  of  the  evaporator  consisted  of  six  differential 
or  state  equations  and  about  sixty  algebraic  equations.  This  latter 
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group  consisted  of  connection  equations  between  the  effects,  property 
correlations  and  heat  transfer  equations.   To  overcome  the  uncertainty 
in  the  empirical  relationships  for  the  inside  film  coefficient  two 
unknown  parameters  were  introduced,  one  for  each  effect.  These  para- 
meters were  estimated  by  correlating  model  predictions  with  data 
collected  on  experimental  runs.  A  nonlinear  least-squares  technique 
was  utilized  to  get  the  best  fit. 

The  algorithm  used  for  the  theoretical  development  of  a  minimum 
time  policy  is  one  in  which  Hamiltonian  minimizations  result  in  control 
policy  updates  on  successive  iterations.  Control  variable  constraints 
and  point  constraints  are  accounted  for  along  the  trajectory.  The 
utility  of  the  algorithm  is  enhanced  by  assuming  a  control  scenario  and 
determining  whether  it  is  optimal  when  compared  to  other  likely  scenarios, 
This  approach  keeps  the  number  of  active  state  and  adjoint  variables  to 
a  minimum  at  any  particular  time  resulting  in  simple  Hamiltonians  and 
less  computational  expense  for  the  integration  of  the  state  and  adjoint 
equations  and  Hamiltonian  minimizations. 

The  minimum  time  algorithm  was  used  on  the  evaporator  model  to 
determine  the  optimal  policy  under  three  different  sets  of  conditions. 
In  the  first  case  it  was  assumed  that  there  was  a  constraint  on  the 
maximum  value  of  the  second  effect  hold-up.  The  second  case  dealt  with 
a  different  set  of  control  variables  in  that  the  feed  rate  to  the  second 
effect  was  assumed  to  be  fixed.  In  the  third  case  the  assumption  of  a 
constrained  maximum  second  effect  hold-up  was  done  away  with.  The 
simulation  results  indicated  that  the  third  case  resulted  in  the  smallest 
start-up  time  with  the  optimal  policy  calling  for  an  overfilling  of  the 
second  effect  followed  by  a  gradual  decrease  in  the  second  effect  hold-up 
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to  the  desired  value  which  took  place  when  boiling  just  started  in 
the  second  effect. 

For  all  three  cases  it  was  found  that  the  control  policy  is 
bang-bang  in  nature  and  that  the  control  switches  occur  at  times  at 
which  the  point  constraints  are  met  on  the  assumed  scenario.  Because 
of  this  the  switching  times  can  be  related  to  the  state  variables 
and  a  feedback  control  policy  is  obtained. 

Experiments  were  run  to  try  out  the  optimal  control  policy  and 
to  test  the  model.  On  an  average,  the  simulations  resulted  in  final 
times  which  were  between  ten  and  fifteen  percent  within  those  obtained 
experimentally.  This  accuracy  was  reasonable  considering  the  experi- 
mental problems  associated  with  hold-up  measurements  and  analog  control 
of  the  hold-ups  and  the  theoretical  problems  associated  with  the  as- 
sumption of  the  heat  transfer  mechanisms. 
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CHAPTER  I 
INTRODUCTION 

Optimal  control  theory  has  been  developed  to  a  fairly  sophis- 
ticated level  in  the  field  of  electrical  engineering.  However,  the 
uses  of  the  theory  and  possible  applications  in  chemical  engineering 
have  been  virtually  unexplored.  The  reasons  for  this  rather  limited 
progress  on  both  the  theoretical  and  applied  fronts  are  many  (Foss, 
1973). 

The  starting  point  in  the  applications  of  control  theory  is 
a  good  dynamic  model  of  the  process.  Most  chemical  processes  have 
been  modeled  poorly  due  to  an  incomplete  understanding  of  the  complex 
interactions  among  numerous  variables.  High  dimensionality  and  non- 
linearities  in  behavior  require  the  use  of  sophisticated  numerical 
techniques  for  the  simulation  and  design  of  control  schemes.  Many 
chemical  processes  have  inherently  large  time  constants  which  make 
them  unsuitable  for  control.  It  is  generally  not  possible  to  make 
all  the  measurements  that  are  required  for  a  feedback  control  scheme, 
and,  in  addition,  measurements  are  subject  to  noise  which  is  not 
readily  filtered  in  nonlinear  systems.  Since  most  of  the  states  are 
not  measurable,  there  is  a  need  for  building  reduced  order  observers 
which  again  is  a  formidable  problem  for  nonlinear  systems. 

In  spite  of  the  above  drawbacks,  quite  a  few  articles  on  the 
subject  of  optimal  control  have  appeared  in  the  chemical  engineering 


literature  during  the  last  five  years.  The  early  investigations  dealt 
with  the  control  of  simplified  lumped  parameter  linear  processes.  The 
linear  quadratic  loss  problem  resulted  in  feedback  control  which  was 
particularly  useful  in  regulatory  control;  i.e.  control  in  the  face  of 
disturbances  (Nieman  and  Fisher,  1970;  Newell  and  Fisher,  1971;  Newell 
et  al.,  1972).  The  arbitrary  nature  of  the  correlation  of  the  weighting 
matrices  to  the  actual  physical  problem  often  makes  the  linear  quadratic 
loss  problem  unrealistic.  Later  investigators  extended  these  techniques 
to  the  control  of  nonlinear  lumped  parameter  systems  by  using  various 
forms  of  linearization  on  the  system  equations  (Weber  and  Lapidus,  1971a, 
1971b;  Siebenthal  and  Aris,  1964;  Tsang  and  Luus ,  1973).  Others  worked 
on  nonlinear  systems  with  one  or  two  control  and  state  variables  (Joffe 
and  Sargent,  1971;  Jackson,  1966).  The  control  of  distributed  parameter 
systems  is  still  in  its  infancy.  Some  investigators  have  reported  sub- 
optimal  control  of  distributed  systems  in  which  some  other  criterion, 
such  as  minimization  of  a  Lyapunov  functional,  is  used  (Vermeychuk  and 
Lapidus,  1973a,  1973b;  Chant  and  Luus,  1968).  Simulated  start-up  studies 
have  been  made  on  plate  distillation  columns  (Pollard  and  Sargent,  1966) 
and  autothermic  reaction  systems  (Jackson,  1966). 

This  work  was  directed  towards  applying  existing  control  theory 
to  a  useful  chemical  process.  The  aim  was  to  study  the  dynamics  and 
control  of  a  simple,  yet  reasonably  complex,  piece  of  equipment  commonly 
found  in  the  chemical  industry.  A  double  effect  evaporator  was  chosen 
as  the  subject  of  study  for  the  above  reason  and  also  because  a  laboratory 
scale  double  effect  evaporator  was  available  for  experimental  work.  The 
study  involved, 
a)  Developing  a  nonlinear  model  for  the  evaporator. 


b)  Estimating  model  parameters  to  fit  experimental  data. 

c)  Developing  a  minimum  start-up  time  control  policy  taking  into 
account  constraints  on  the  state  and  control  variables  and  putting 
the  optimal  policy  in  feedback  form  in  terms  of  switching  times. 

d)  Experimentally  determining  the  effect  of  the  policy. 

This  approach  differs  from  previous  ones  in  a  few  respects. 
The  model  is  highly  nonlinear  and  is  treated  as  such.  No  linearization 
is  resorted  to  as  start-up  involves  large  changes  in  the  state  variables 
and  linearized  equations  would  be  inaccurate.  The  mathematics  involved 
in  obtaining  the  minimum  time  policy  is  simplified  as  the  approach 
adopted  presupposes  a  start-up  scenario  and  then  verifies  that  it  is 
optimal.  The  algorithm  leading  to  the  optimal  policy  handles  con- 
straints on  control  and  state  variables  in  a  logical  fashion  by  directly 
holding  the  state  or  control  variable  on  the  constraint  and  changing 
the  equation  set  and  its  solution  procedure  as  a  result.  This  avoids 
the  use  of  penalty  function  methods  and  the  like.  Finally,  the  control 
policy  is  experimentally  verified. 

The  minimum  time  objective  was  chosen  primarily  from  the  point 
of  view  of  economics.  Control  costs  during  start-up  are  minimal  com- 
pared to  the  start-up  time  in  batch  processes.  Reducing  the  start-up 
time  results  in  reduced  down  time  thus  improving  cycle  efficiency  and 
increasing  profits.  The  food  industry  is  an  example  of  an  industry 
which  must  shut  down  frequently  to  have  the  processing  equipment  cleaned. 
Orange  juice  is  concentrated  in  multiple  effect  evaporator  systems,  and 
these  systems  are  cleaned  about  three  times  a  day.  A  second  reason  for 
minimum  time  start-up  was  more  specific  to  the  ultimate  use  of  the  parti- 
cular double  effect  evaporator  investigated.  It  is  to  be  used  in  an 


undergraduate  laboratory  experiment  in  computer  control,  and  past 
experiences  indicated  that  it  took  a  very   long  time  to  bring  it  to 
steady  state  under  manual  control.  Thus,  in  order  to  reduce  the 
start-up  time  and  consequently  to  reduce  the  amount  of  on-line  com- 
puter time  for  steady  state  observations,  it  was  imperative  to  have 
start-up  in  a  minimum  time. 

Chapter  II  contains  a  description  of  the  experimental  evapo- 
rator, the  instrumentation  and  the  interfacing  equipment  with  the 
IBM  370/165  (which  is  the  main  computer  on  campus).  Chapter  III  deals 
with  the  building  of  a  dynamic  model  for  the  evaporator  and  also  the 
estimation  of  parameters  from  experimental  data.  A  derivation  of  the 
optimal  control  algorithm  is  given  in  Chapter  IV.  It  also  contains  the 
simulated  and  experimental  results  of  the  application  of  the  control 
algorithm.  Some  comments  and  proposals  for  further  work  are  given  in 
Chapter  V.  Appendix  A  contains  all  the  heat  transfer  equations  which 
supplement  the  main  model  equations  in  Chapter  III.  A  listing  and 
description  of  the  computer  program  implementing  the  optimal  control 
algorithm  is  the  subject  of  Appendix  B. 


CHAPTER  II 
DESCRIPTION  OF  EVAPORATOR  AND  COMPUTER  INTERFACING  EQUIPMENT 

II.  1  Evaporator  Layout  and  Description 

The  double  effect  evaporator  is  located  in  the  unit  operations 
laboratory  of  the  chemical  engineering  department.  Figure  2.1  is  a 
schematic  of  the  double  effect,  showing  the  arrangement  of  the  two 
effects,  EV1  and  EV2,  and  the  basic  process  and  vapor  lines.  Note 
that  backward  feed  is  used;  that  is,  the  vapor  flow  and  process  fluid 
flow  are  in  opposite  directions. 

The  first  effect  is  a  long  tube  vertical  (LTV)  evaporator. 
It  contains  3  tubes,  each  9  feet,  6  inches  (2.90  m)  long  and  1  inch 
(0.0254  m)  O.D.  with  heating  steam  at  about  20  psig  (2.39  bar)  on  the 
outside  of  the  tubes.  The  process  fluid  flows  upward  through  the 
tubes  either  by  natural  or  forced  circulation.  The  latter  method  is 
almost  always  used  because  of  the  increased  heat  transfer  coefficients 
obtainable.  The  pressure  on  the  process  side  is  at  or  slightly  above 
atmospheric. 

The  mixture  of  process  fluid  and  vapor  formed  in  the  first 
effect  enters  a  vapor-liguid  separator,  SE1 ,  which  is  at  the  same 
pressure  as  the  first  effect.  The  liquid  is  drawn  off  the  bottom  of 
the  separator  and  is  recirculated  back  into  the  first  effect  by  pump 
PU3  after  some  liquid  product  is  withdrawn.  Fresh  feed  to  the  first 


K  A  ,\  A  A  A  A 

V  V         m  /  v  V  V 


UJ   Q 
Q 

o  1— 


8 


effect  is  pumped  by  pump  PU2  from  the  second  effect.  The  vapor  from 
separator  SE1  is  used  as  the  steam  input  in  the  second  effect.  This 
leads  to  steam  economy  as  one  pound  of  heating  steam  used  in  the  first 
effect  should  evaporate  more  than  one  pound  of  water  from  the  first 
and  second  effects  combined. 

The  second  effect,  EV2,  is  a  calandria  type  effect  in  which 
there  are  15  tubes,  each  2  feet,  4  inches  (0.711  m)  long  and  1  inch 
(0.0254  m)  O.D.  The  effect  also  has  a  2  inch  (0.0508  m)  O.D.  central 
downtake.  Heat  transfer  is  by  natural  convection  only,  resulting  in 
much  lower  heat  transfer  coefficients  compared  to  the  first  effect. 
Fresh  preheated  feed  is  pumped  into  the  bottom  of  the  second  effect 
from  the  feed  tank  by  pump  PU1 .  The  heating  medium  is  the  vapor  from 
the  first  effect  on  the  outside  of  the  tubes.  Above  the  calandria  is 
a  vapor  body  which  separates  the  vapor  from  the  liquid.  The  vapor  is 
drawn  into  a  condenser,  CD1  ,  by  means  of  a  vacuum  produced  by  a  steam 
jet  ejector.  The  ejector  maintains  the  pressure  on  the  process  side 
in  the  effect  at  around  10  inches  mercury  vacuum  (0.675  bar). 

The  vapor  condensate  from  the  first  effect  is  collected  in  tank 
Tl  and  that  from  the  second  effect  is  collected  in  tank  T2,  both  of 
which  are  maintained  at  a  vacuum  by  the  same  steam  jet  ejector. 

1 1. 2  Operating  Notes 

There  were  a  few  precautions  which  had  to  be  observed  during 
operation. 

1)  The  feed  rate  was  kept  at  around  2-3  gpm  (0.12  to  0.18  kg/s). 

2)  The  recirculation  rate  in  the  first  effect  was  kept  at  a  maximum 
of  15  gpm  (0.9  kg/s).  A  higher  rate  caused  entrainment  of  liquid  with 


the  vapor  in  the  vapor-liquid  separator  SE1 .  This  separator  has  no 
baffling  of  any  kind  and  is  very   inefficient  at  high  flow  rates. 

3)  To  avoid  cavitation  in  the  recirculation  pump  PU3,  care  was  taken 
to  see  that  the  vertical  suction  leg  from  the  separator  to  the  re- 
circulation pump  was  always  filled  with  liquid.  This  was  particularly 
critical  when  the  pump  was  first  started.  Incomplete  filling  of  the 
vertical  leg  led  to  pulsatinq  flows  resulting  in  large  upsets  in  the 
evaporator  operation.  A  recirculation  rate  higher  than  15  gpm  (0.9 
kg/s)  also  caused  a  high  discharge  head  on  pump  PU3,  much  higher  than 
the  maximum  discharge  head  on  pump  PU2  (which  is  of  a  smaller  capacity), 
eliminating  all  flow  of  fresh  feed  to  the  first  effect. 

4)  The  liquid  level  in  the  second  effect  was  maintained  around  the 
top  of  the  tubes  for  best  utilization  of  the  heat  transfer  area. 

II. 3  Evaporator  Instrumentation 

As  part  of  this  work  the  evaporator  piping  had  to  be  modified 
to  accomodate  the  instrumentation  required  for  control.  The  work  con- 
sisted mainly  of  installing  pneumatic  control  valves,  orifices,  pressure 
taps,  thermocouples  and  extra  manual  valves.  Figures  2.2  and  2.3  show 
the  detailed  instrumentation  of  the  evaporator.  The  legend  of  Fiqure 
2.2  also  applies  to  Figure  2.3. 

Three  pneumatic  control  valves,  CV1 ,  CV2  and  CV3,  were  installed 
in  the  feed,  inter-effect  and  recirculation  lines  respectively.  CV2 
and  CV3  were  normally  closed  (air-to-open)  valves  and  were  installed  in 
bypasses  on  the  lines,  whereas  CV1  was  a  normally  open  valve  and  was  in- 
stalled in  the  feed  line  as  such.  The  purpose  of  the  by-passes  was  to 
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allow  for  complete  or  partial  manual  control  of  experimental  runs 
when  desired. 

Flow  rates  were  measured  with  square-edged  orifices,  0R1 , 
0R2,  0R3  and  0R4,  installed  in  the  feed,  inter-effect,  recirculation 
and  product  lines  respectively.  The  pressure  drop  across  an  orifice 
is  indicative  of  the  flow  through  it. 

Liquid  levels  (proportional  to  hold-ups)  were  measured  by 
taking  the  difference  in  total  pressure  between  the  bottom  and  top 
of  each  of  the  two  effects.  Pressure  taps  were  installed  at  both 
ends  of  the  sight  glasses  for  this  purpose.  The  upper  taps  were 
also  used  to  measure  absolute  pressure  in  the  effects. 

Temperatures  were  measured  by  jacketed  copper-constantan 
thermocouples,  TCI,  TC2,  TC3  and  TC4.  These  were  installed  in  the 
feed  line,  at  the  exit  of  the  second  effect,  at  the  entrance  to  the 
first  effect,  and  in  the  steam  chest  of  the  first  effect,  severally. 

All  liquid  lines  from  the  pressure  taps  and  air  lines  to  the 
valves  were  brought  to  a  central  panel  in  the  front  of  the  evaporator 
with  1/4-inch  poly-flo  tubing.  Quick-connect  fittings  were  used  at 
the  panel  so  that  leads  to  the  interfacing  equipment  could  be  con- 
nected quickly  when  required.  The  thermocouple  wires  also  terminated 
with  special  thermocouple  outlets  at  the  panel. 

II. 4  Transducing  and  Controlling  Equipment 

The  transducing  and  controlling  equipment  was  installed  in  a 
19-inch  relay  rack  on  casters.  All  air  and  liguid  lines  were  of  poly- 
flo  tubing  with  quick-connect  fittings.  This  made  the  rack  very  ver- 
satile as  it  can  be  moved  to  a  number  of  different  pieces  of  equipment 
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if  desired.  A  layout  of  the  cabinet  is  shown  in  Figure  2.4. 

The  pneumatic  controllers  have  adjustable  proportional  and 
reset  action,  motorized  set  point  control,  and  indication  facilities. 
One  Fischer  and  Porter  model  51  and  three  Taylor  model  662R  control- 
lers were  installed.  The  set  point  motor  of  the  Fischer  and  Porter 
controller  operates  on  a  pulse  train  input  and  Taylor  controllers  on 
a  24-volt  DC  signal.  The  controllers  are  equipped  with  feedback 
potentiometers  which  indicate  their  set  point  positions.  The  pneu- 
matic input  signal  range  to  each  controller  (from  the  DP  cells)  is 
3-15  psig  and  that  of  the  output  pressure  to  the  associated  valve 
is  also  a  3-15  psig  signal . 

The  EMF  to  pneumatic  converters  (not  used  in  the  current  ex- 
periments) are  Foxboro  model  33A  converters.  They  can  transduce 
either  a  millivoltage  or  voltage  signal  into  a  3-15  psig  pneumatic 
signal . 

The  differential  pressure  (DP)  cells  used  are  Foxboro  Model 
13A  DP  cells.  The  adjustable  range  of  the  differential  input  signal 
is  0-500  inches  water  and  the  proportional  pneumatic  output  is  in  the 
3-15  psig  range.  These  DP  cells  were  used  to  transduce  the  pressure 
drops  across  the  orifices  and  the  pressure  differences  corresponding 
to  the  liquid  levels.  The  outputs  were  thus  proportional  to  the  flow 
rates  or  to  the  liquid  levels  (hold-ups).  The  DP  cells  were  also 
used  for  measuring  absolute  pressures  by  venting  the  high  pressure 
side  when  measuring  pressures  above  atmospheric.  The  output  pressure 
in  this  case  was  proportional  to  the  vacuum  or  above-atmospheric 
pressure. 
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II. 5  IBM  1070  Interface 

An  IBM  1070  interface  was  used  which  consisted  of  an  IBM  1071 
central  terminal  unit,  two  IBM  1072  multiplexing  units,  an  IBM  1073 
Model  3  digital  to  pulse  converter  and  an  IBM  1075  decimal  display. 
The  interface  with  all  the  auxiliary  equipment  resides  in  3  relay  racks 
as  shown  in  Figure  2.4.  These  cabinets  were  originally  assembled  by 
R.  C.  Eschenbacher  and  are  described  in  his  Ph.D.  thesis  (Eschenbacher , 
1970).  However,  as  part  of  this  work,  some  of  the  equipment  had  to  be 
rewired  to  accomodate  pulse  duration  outputs  and  most  of  the  relays 
had  to  be  rewired  and  chanqed  to  double  coil  relays  in  order  to  isolate 
the  IBM  equipment  from  the  user  circuitry.  The  entire  set  up  is  de- 
scribed in  detail  in  the  GIPSI  (General  .Interface  for  Process  Systems 
Instrumentation)  hardware  manual  (GIPSI,  197  3)  and  is  summarized  very 
briefly  here. 

1)  Input 

The  pressure-voltage  transducers  converted  the  3-15  psig  air 
signals  from  the  DP  cells  to  0-5  volt  DC  signals  which  were  fed  as 
analog  inputs  to  the  1070.  The  millivoltage  thermocouple  signals 
were  also  fed  as  analog  inputs  through  the  special  thermocouple  in- 
put feature  in  the  1070.  The  unit  also  has  a  facility  for  digital 
input  which  was  not  used. 

Another  convenient  input  facility  frequently  used  was  a  form 
of  digital  input  throuqh  pre-designated  demand  functions  which  were 
dialled  into  the  system  through  rotary  switches. 

2)  Output 

Output  from  the  1070  was  in  digital  and  in  pulse  form.  The 
digital  output  was  used  mainly  to  ring  a  bell  to  alert  the  operator 
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to  possible  alarm  conditions  in  the  hardware  and  software.  The 
pulsed  output  was  obtained  from  the  1073  and  was  used  to  move  the  set 
points  on  the  controllers.  The  digital  to  pulse  converter  (1073) 
has  outputs  in  the  form  of  a  pulse  train  as  well  as  a  duration  pulse. 

The  1075  decimal  display  is  a  feature  which  utilizes  digital 
outputs  and  was  used  to  display  particular  variable  values  or  error 
codes. 
3)  Process  Alerts 

The  1070  interface  is  linked  with  the  IBM  370/165  computer  on 
campus.  The  seven  process  alerts  attached  to  a  process  alert  (PA) 
bus  in  the  1070  provide  a  hardware  interrupt  capability  of  the  com- 
puter by  the  process.  The  software  issues  a  conditional  read  of  the 
1070  terminal  to  the  computer.  It  is  then  in  a  hardware  wait  stage. 
When  the  PA  bus  is  activated,  by  one  of  the  process  alerts  on  the 
1070,  the  IBM  370/165  computer  senses  the  closure  and  reactivates 
the  software  which  then  determines  which  PA  was  set.  The  software 
then  resets  the  PA  and  executes  the  program  associated  with  the  PA. 
Process  alert  1  has  the  added  facility  of  being  set  automatically  by 
a  hardware  poller  on  which  the  timing  is  adjustable.  This  enables 
one  to  have  PA1  periodically  and  automatically  set  after  a  predeter- 
mined time  interval  has  elapsed. 

Figure  2.5  is  a  schematic  illustrating  the  flow  of  information 
among  the  various  hardware  components  of  the  experiment. 

1 1. 6  Software 

The  GIPSI  software  package  was  originally  written  by  R.  C. 
Eschenbacher.  Version  2,  the  version  used,  was  written  by  L.  A.  Delgado 
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( G IPS  1 ,  1973)  and  was  marginally  modified  as  part  of  this  work  to 
extend  its  capabilities.  It  is  written  entirely  in  Fortran  (with 
the  exception  of  certain  input/output  routines  in  BTAM  provided  by 
IBM),  has  an  extensive  debugging  facility,  and  an  extensive  error 
handling  facility  to  flag  user  software  and  hardware  errors.  A 
simplified  flow  chart  is  reproduced  from  Westerberg  and  Eschenbacher 
(1971)  and  is  shown  in  Figure  2.6.  It  is  described  in  greater  detail 
in  Westerberg  and  Eschenbacher  (1971)  and  GIPSI  (1973). 

The  heart  of  the  software  is  the  concept  of  the  execute  and 
delay  stacks.  When  a  PA  is  set,  it  is  identified  and  the  program 
(or  programs)  associated  with  it  are  stacked  by  the  program  stacker 
in  the  order  of  priority  on  the  execute  stack.  Control  then  passes 
to  the  Execute  subprogram  which  then  examines  the  execute  stack  and 
passes  control  one  at  a  time  to  the  programs  that  are  due  for  execu- 
tion. If  the  sequence  began  with  PA1 ,  its  response  program  CLOCK 
removes  programs  from  the  delay  stack  if  their  delay  time  has  expired 
and  puts  them  on  the  execute  stack.  The  Execute  subprogram  then  finds 
additional  programs  on  the  execute  stack  which  it  continues  to  remove 
and  cause  to  be  executed.  This  is  done  until  the  execute  stack  is 
empty  whence  control  returns  to  the  PA  handler  which  issues  a  condi- 
tional read  and  the  IBM  370/165  again  waits  for  a  PA  to  be  set  to 
start  the  cycle  again.  All  delay  times  are  compared  to  the  computer 
clock.  Data  on  program  priorities,  delay  time,  etc.  are  specified 
in  the  Program  Descriptive  Data. 

All  the  user  has  to  do  is  to  provide  his  specific  programs 
associated  with  the  various  process  alerts,  the  program  descriptive 
data  for  all  his  programs,  and  a  subroutine  GOTO  which  the  Execute 
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subprogram  uses  to  pass  control  to  CLOCK  and  the  user  and  system 
subprograms. 

The  computer  costs  are  extremely  low  when  based  mainly  on 
central  processing  unit  (CPU)  time.  The  software  utilizes  very 
little  CPU  time.  Typical  costs  are  in  the  range  of  $3-5  an  hour 
provided  no  elaborate  computations  are  called  for  in  the  user  pro- 
grams. However,  costs  for  core  residency  charges  dominate  as  the 
basic  software  package  requires  around  20,000  words  or  80,000 
bytes  of  core. 


CHAPTER  III 
DYNAMIC  MODEL  AND  PARAMETER  ESTIMATION 

The  dynamic  modeling  of  multiple-effect  evaporators  has  been 
extensively  investigated  in  recent  years  at  the  University  of  Alberta 
(Andre  and  Ritter,  1968),  (Newell,  1970).  In  simulation  and  experi- 
mental work  high  order,  linear  models  have  been  found  to  be  satisfactory. 
However,  linear  models  are  not  realistic  when  the  operating  conditions 
change  drastically  as  in  start-up.  In  the  first  part  of  this  chapter 
a  nonlinear,  first  order,  lumped  parameter  model  is  proposed.  The 
first  order  and  lumped  parameter  nature  of  the  model  was  resorted  to 
for  two  main  reasons: 

1)  The  model  was  simple  and  adequately  described  the  data. 

2)  The  model  was  used  to  devise  an  optimal  control  policy  for  minimum 
time  start-up.  Optimal  control  theory  has  been  rigorously  developed 
for  lumped  parameter  systems  and  its  extension  to  distributed  systems 
has  not  yet  been  extensively  investigated. 

In  addition,  the  model  presented  here  takes  into  account  heat 
transfer  dynamics  from  the  viewpoint  of  film  coefficients.  Although 
this  leads  to  complicated  algebraic  equations,  it  has  the  advantage  of 
leading  to  a  better  understanding  of  the  heat  transfer  dynamics.  It 
also  gives  rise  to  two  constant  correction  parameters.  The  necessity 
for  these  parameters  is  due  to  the  uncertain  coefficients  that  are  used 
in  the  film  coefficient  equations.  The  second  part  of  this  chapter 
deals  with  the  estimation  of  these  parameters  to  fit  the  experimental  data 
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HI  .1  Dynamic  Model 

The  dynamic  model  is  a  collection  of  the  material  and  energy 
balances  for  each  effect.  For  a  double  effect  evaporator  concen- 
trating a  solution  with  one  major  solute,  there  are  two  material 
balances  (one  for  the  solution  and  one  for  the  solute)  and  one  energy 
balance  for  each  effect,  giving  rise  to  a  total  of  six  dynamic  or 
state  equations  for  the  two  effects.  In  addition,  there  are  dynamic 
equations  for  the  vapor  phases  and  metal  but  the  time  constants  of 
these  are  negligible  compared  to  the  six  mentioned  earlier  (Andre, 
1968)  so  that  these  dynamic  equations  could  be  reduced  to  be  alge- 
braic equations.  This  procedure  of  setting  the  derivatives  of  the 
equations  with  small  time  constants  to  zero  reduces  the  order  of 
the  system.  The  full  model  will  be  presented  here.  In  later  chapters, 
appropriate  simplifications  will  be  applied  as  some  of  the  model 
states  are  held  fixed  (for  example,  as  boiling  does  or  does  not  take 
place).  A  summary  of  all  the  assumptions  made  is  presented  at  the 
end  of  the  model.  Refer  to  Figure  3.1  for  the  symbols  used  for  the 
flows,  hold-ups  and  temperatures.  There  is  also  a  foldout  nomenclature 
list  on  page  148. 
III. 1.1  State  Equations 

1)  First  and  second  effect  hold-ups,  H-,  and  bL. 

dT-=W12-V21  -«11  -W01  (3J) 

dl!0 

3T  '  WF  "  V02  "  W12  (3'2) 

2)  First  and  second  effect  enthalpies. 
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Since  the  evaporator  was  to  be  used  ultimately  for  the  con- 
centration of  dilute  solutions,  it  was  assumed  that  there  would  be 
no  boiling  point  elevations  in  the  effects.  Also,  perfect  mixing 
is  assumed  which  would  be  close  to  the  case  for  small  hold-ups  and 
dilute  solutions. 
d(H1h  ) 

-at— =  wi2hi2  -  v2ihI-  <wn  +1Whi  +ni 

where  Q,  is  the  heat  transferred  from  the  steam.  Simplifying  this 
using  equation  (3.1),  we  get 

dh1   , 

dT"=  M7[H12(h12  "  V  +  V21<hl  "  hl>  +  V         (3-3 

Similarly,  an  energy  balance  on  the  second  effect  under  the 

same  assumptions  gives  rise  to 

dh 

3)  First  and  second  effect  solute  material  balances. 

Again,  assuming  perfect  mixing  we  have  for  the  first  effect 
solute, 


ar =  ft-  EwF(hF  -  h2)  +  v02(h2  -  hp  +  o2]  (3.4) 


d(H1Cj 

-dT~  =W12C12-   <wn   +  l'WCl 


Simplifying  this  with  equation  (3.1)  we  have 

dCl        1 

dT~-  TTj"  CW12^C12   "  Cl}   +  V21C1]  (3.5) 

Similarly,  a  balance  on  the  second  effect  solute  yields, 
dC 

ar-TfC^F  -  c2^ +  vo2y  (3.6) 
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III. 1.2  Connection  Equations 

In  addition  to  the  state  equations  listed  above,  there  are 
algebraic  equations  which  arise  due  to  the  mixing  of  the  two  streams 
between  the  second  and  first  effects.  One  energy  and  two  material 
balances  describe  the  mixing  as  follows: 

W12  =  1^2  +  Wn  <3'7) 

w12c12  =  w-i2c2  +  lJnci  (3,8) 

W]2h12  =  W]2h2  +  l!11h1  (3.9) 

III. 1.3  Heat  Transfer  Equations 

These  equations  arise  in  computing  the  terms  0]  and  0.2  that 
arise  in  the  enthalpy  equations  (3.3)  and  (3.4). 

The  heat  transfer  rate  Q,  is  a  function  of  the  steam  temper- 
ature, the  first  effect  temperature,  the  inside  and  outside  film 
coefficients,  the  wall  resistance  including  fouling  and  the  heat 
transfer  area  in  the  first  effect.  The  inside  film  coefficient  is 
a  function  of  the  flow  rate  through  the  tubes,  the  vaporization, 
the  inner  wall  and  bulk  temperatures  and  the  entrance  temperature. 
The  heat  transfer  mechanism  is  initially  simple— a  combination  of 
the  Dittus-Boelter  equation  for  the  inside  and  the  Nusselt  equation 
for  the  condensing  steam.  However,  when  boiling  takes  place  two- 
phase  heat  transfer  occurs  because  of  the  vapor  formed.  The  complete 
boiling  mechanism  is  a  topic  for  further  investigation.  Approximate 
correlations  were  obtained  from  (Fair,  1960,  1963a,  1963b),  (Hughmark 
1969)  and  more  recently  from  (Tong,  1965).  The  complete  list  of 
equations  leading  to  the  determination  of  C^  from  the  state  variables 
and  flow  rates  is  given  in  Appendix  A.  Due  to  the  uncertainty  in 
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the  empirical  equations  which  predict  the  inside  film  coefficients, 
a  parameter  6-.  was  introduced  in  the  overall  heat  transfer  equations 
(A. 17)  and  (A. 40).  It  is  assumed  that  the  outside  film  coefficient 
is  predicted  by  the  Nusselt  equation,  (A. 14)  and  (A. 28)  to  a  rea- 
sonably high  degree  of  accuracy  as  is  borne  out  later  by  experiment. 
It  is  also  assumed  that  the  parameter  9-,  has  two  different  values 
depending  upon  the  heat  transfer  mechanism  in  the  first  effect.  This 
depends  upon  the  stage  of  start-up  as  follows: 

1)  9-,  =  9-,  ,  when  the  first  effect  liquid  is  being  heated.  The 

i  a 

dominant  equation  for  the  inside  film  coefficient  is  solely  the  Dittus- 
Boelter  equation  (A. 16). 

2)  9,  =  9-,,,  when  the  liquid  in  the  first  effect  is  boiling.  The 
inside  film  coefficient  is  a  combination  of  many  factors  including  a 
coefficient  due  to  nucleate  boiling  (A. 38)  and  a  two-phase  convective 
coefficient  (A. 33). 

The  equation  for  the  first  effect  heat  transfer  rate  is  given 

in  functional  form  as: 

Q1  =  Q1(Ts,T1,T12,HTW12,V21,e1)  (3.10) 

where  the  temperatures,  T,  are  functions  of  the  enthalpies,  h,  in  the 

form 

T  =  f(h) 

Equation  (3.10)  has  implicitly  used  the  fact  that  the  overall 
coefficient,  film  coefficients  and  heat  transfer  area  are  functions  of 
temperatures  (enthalpies)  and  hold-ups. 

Strictly  speaking,  the  equations  describing  the  steam  (vapor) 
temperatures  or  enthalpies,  h  ,  h,  and  h,>  are  differential  equations, 
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dp, 
Vi  HT  =  V21  "  Wc2 

Vn  -^-L-=  V21hT  -  Wc2hc2  -  ^2 

where  V]  is  the  volume  of  the  vapor  space  in  the  first  effect,  the 
vapor-liquid  separator  and  the  tubes  of  the  second  effect.  Note  that 
two  assumptions  have  been  made  here-the  steam  (vapor)  is  saturated 
and  there  is  no  subcooling  of  the  condensate. 

However,  it  has  been  shown  by  Andre  and  Ritter  (1968)  that 
the  response  rate  of  the  steam  enthalpy  is  negligible  compared  to 
that  of  the  hold-up,  concentration  and  enthalpy  equations  (3.1)  to 
(3.6)  in  the  two  effects.  The  differential  equations  describing 
the  steam  density  and  temperature  are  so  replaced  with  the  steady 
state  equations 

V21  =  Wc2 
^  Q2  =  V2lK  -hc2> 

°r  Q2=V21X1  (3>11) 

where  ^  =  f{J^ 

The  heat  transfer  rate  in  the  second  effect,  Q2»  is  also  a 
function  of  the  film  coefficients,  wall  resistance  including  fouling, 
area  and  temperatures  in  the  second  effect.  The  heat  transfer 
mechanism  is  purely  natural  convection.  Here  again,  it  is  assumed 
that  the  Musselt  equation  (A. 47)  is  reasonably  accurate  in  predicting 
the  outside  film  coefficient.  The  inside  film  coefficient  is  pre- 
dicted by  the  natural  convection  equation  (A. 51)  and  the  overall 
coefficient  (A. 52)  has  an  undetermined  parameter  62  which  again  can 
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have  two  values.  One  value  (0„  )  is  for  the  heating  of  the  liquid 
and  the  other  (8ok)  is  f°r  the  boiling  of  the  liquid  in  the  second 
effect.  The  functional  form  for  Q?  is: 

Q2  =  Q2(TrT2,TF,H2,e2)  (3.12) 

where  the  temperatures  have  been  determined  from  the  corresponding 
enthalpies. 

Note  that  there  are  no  liquid  flow  and  vapor  flow  terms  here 
as  the  natural  convection  overall  coefficient  is  not  a  function  of 
these  variables . 

A  subroutine  called  HEAT  has  been  written  to  calculate  the 
heat  transfer  rates  Q-,  and  Q?  from  the  temperatures,  hold-ups  and 
flow  rates.  It  is  included  in  Appendix  B.  The  rates  0-,  and  Q„  are 
found  by  using  all  the  heat  transfer  equations  in  Appendix  A.  The 
inner  and  outer  wall  temperatures  that  figure  in  the  film  coefficient 
calculations  are  unknown.  These  temperatures  are  initially  guessed 
and  the  film  coefficients  are  calculated.  The  wall  temperatures  are 
adjusted  until  the  equations  predict  the  same  heat  transfer  rate 
per  unit  area  across  the  inside  and  outside  films  and  the  wall.  With 
the  final  wall  temperatures,  the  film  coefficients  and  heat  transfer 
rates  are  estimated  using  the  appropriate  equations  depending  upon 
the  nature  of  boiling  and  the  mechanism. 
III. 1.4  Decision  Variables 

Equations  (3.1)  to  (3.6)  are  the  differential  equations  and 
(3.7)  to  (3.12)  are  the  algebraic  equations  describing  the  dynamic 
response  of  the  evaporator.  After  enumerating  the  number  of  variables 
and  the  number  of  equations  it  is  found  that  there  are  eight  more 
variables  than  there  are   equations.  These  eight  decision  variables 
are  chosen  in  a  natural  way  making  them  manipulative  or  control 
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variables.  These  are: 

1 )  feed  flow  rate,  Wp 

2)  feed  temperature,  Tp 

3)  feed  concentration,  Cp 

4)  inter-effect  flow  rate,  W^2 

5)  recirculation  in  first  effect,  W-j -| 

6)  product  flow  rate,  Wq-| 

7)  steam  temperature  in  first  effect,  Ts 

8)  total  pressure  in  second  effect,  P2 
III. 1.5  Assumptions 

A  summary  of  the  assumptions  made  in  writing  the  model  is  pre- 
sented here. 

1)  Time  responses  of  the  vapor  phase  and  the  tube  walls  are  negligible 
compared  to  that  of  the  liquid  phase.  This  results  in  simpler  algebraic 
equations  for  the  vapor  phase  and  tube  walls  and  also  decreases  the 
dimensionality  of  the  model. 

2)  The  vapor  is  saturated  and  is  in  equilibrium  with  the  liquid  at  the 
same  temperature. 

3)  Condensate  on  the  vapor  side  of  the  first  and  second  effects  is  not 
subcooled  and  condensate  hold-up  is  negligible. 

4)  Boiling  point  elevations  due  to  the  presence  of  solute  in  the  two 
effects  is  negligible.  This  is  justified  in  the  case  of  dilute  solutions 

5)  There  is  perfect  mixing  in  the  two  effects  resulting  in  lumped 
parameter  concentration  and  heat  transfer  equations. 

6)  The  heat  transfer  mechanism  in  the  first  effect  is  single  phase  con- 
vection followed  by  two-phase  convection  and  nucleate  boiling  and  that 
in  the  second  effect  is  natural  convection. 
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7)  Slug  flow  is  the  predominant  flow  pattern  in  the  first  effect  when 
boiling  takes  place. 

8)  Heat  losses  are  negligible. 

9)  The  inaccuracy  of  the  heat  transfer  rates  is  due  mainly  to  the 
uncertainty  of  the  inside  film  coefficient  leading  to  undetermined 
parameters  to  correct  for  the  inside  coefficients  alone. 

III. 2  Parameter  Estimation 

II 1 .2.1  Stochastic  versus  Deterministic  Estimation 

An  extensive  review  of  parameter  estimation  technigues  in 
differential  eguations  is  available  in  Nieman  et  al . ,  (1971).  In 
the  deterministic  case  the  simplest  and  most  effective  method  is  a 
least  sguares  fit.  The  problem  is  stated  as  follows: 
Given  the  state  eguations 

X  =  f[X(t),  U(t),  e(t)] 
where  U(t)  are  the  control  or  manipulative  variables  and  o(t)  are   the 
parameters  to  be  estimated  from  experimental  data  Y(t).  The  obser- 
vations are  related  to  the  states  and  controls  by 

Y(t)  =  h[X(t),  U(t)] 
The  problem  is  to  determine  the  parameters  o(t)  such  that  the  model 
"best  fits"  the  given  experimental  data  Y(t). 

Assuming  that  the  parameters  are  constant,  as  in  the  case  of 
the  evaporator,  O(t)  =  0,  the  problem  can  be  reduced  to  a  least-sguares 

estimation 

N       .  2 
Min  I     (Y.  -  Y.) 
G   i=l 


where  Y-  is  a  calculated  observation  at  time  t-  (i=l,...,N)  and  Y.  is 
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the  actual  observation. 

The  alternative  to  the  above  deterministic  estimation  is  the 
problem  of  stochastic  estimation.  It  seemed  that  the  model  would  fit 
the  data  far  better  if  the  parameters  0  were  updated  as  each  measure- 
ment was  made.  Further,  if  the  states  and  measurements  were  subject 
to  process  and  measurement  noise  it  would  be  necessary  to  estimate  the 
states  using  a  nonlinear  or  a  linearized  Kalman  filter.  It  was  apparent 
all  along  that  this  would  require  an  appreciable  amount  of  computer 
time;  however,  the  estimation  technique  was  investigated  off-line  from 
an  academic  viewpoint. 

The  approach  follows  that  of  Padmanabhan  (1970).  The  model 
is  assumed  to  have  process  and  measurement  noise  v(t)  and  w(t) 
X(t)  =  f(X(t),  t)  +  V(t) 
Y(t)  =  h(X(t),  t)  +  W(t) 
where  V(t)  and  W(t)  are  white  Gaussian  noise  sequences  with  zero  mean 
and  covariances  Q(t)  and  R(t)  respectively.  Mote  that  the  control 
vector  U(t)  is  expressed  in  terms  of  the  state  vector  X(t)  in  the  state 
and  observation  equations  above.  To  the  state  equations  could  be 
augmented  the  parameter  equations 

O(t)  -  0  +  V (t) 
thus  treating  the  parameters  0  as  states. 

The  problem  reduces  to  estimating  the  states  X(t)  and  the 
parameters  O(t)  from  the  experimental  data  Y(t).  A  recursive  estimate 
X(t,/t, )  representing  an  estimate  of  X  at  time  t^  based  on  all  the 
data  collected  until  time  tk  is  given  by 

$(tk/tk)  =f(x(tk/tk),  tk)  +  P(tk)z(tk/tk) 
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with  X(0/0)  =  \j  =  initial   estimate  of  X(t  ) 

where  Z(ytk)  =  h^R(tk)_1  (Y(tk)   -  h(X(ytk),  tR)) 

and  the  covariance  of  the  estimate  P(t. )  satisfies  a  matrix  Ricatti 
equation 

HP         ^  ^T         ^ 

&  (tk)  ■  f/(tk)  +  P(tk)fJ  +  P(tk)ZxP(tk)  ♦  Q(tk) 


with  P(0)  =  tt  =   initial  covariance  of  estimate  X(tQ).  In  the  linear 
case  the  covariance  equation  can  be  integrated  off-line  and  stored  for 
use  by  the  state  estimation  equations. 

To  study  the  effectiveness  of  the  scheme  an  example  was  run 
off-line  in  which  it  was  assumed  that  the  first  and  second  effect  hold- 
ups and  the  first  effect  temperature  were  constant.  It  was  desired  to 
estimate  the  second  effect  temperature  and  the  parameter  e1  in  the 
equation 

dh?   1 

^.  =  ^[wF(hF-h2)  +  e'Q2] 

while  measurements  were  obtained  on  the  second  effect  temperature 
T2  =  f(h2) 
The  state  equation  and  the  Ricatti  equations  were  integrated 
for  eight  minutes  of  real  time  and  this  took  100  seconds  of  IBM  370/165 
CPU  time.  The  results  for  the  first  four  minutes  were  as  follows: 


Time 

x4  .measured 
118.36 

X4 

x4,  predicted 

e' 

30.7 

1.0 

31.7 

122.63 

121.4 

118.73 

0.08 

32.7 

124.32 

124.1 

122.37 

0.18 

33.7 

125.82 

126.76 

125.05 

0.19 

34.7 

126.62 

127.34 

126.35 

0.17 
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This  example  showed  that  it  was  not  practical  to  try  on-line 

stochastic  estimation  for  a  problem  of  this  nature  especially  when 

there  are  more  state  variables.  Because  of  this  the  deterministic 
least-squares  estimate  was  resorted  to  and  the  results  obtained  were 
acceptably  good. 

To  account  for  noisy  flow  variables  a  linear  filter  was  used 
on  the  flow  measurements  in  the  form 

^(j)  =  au.(j  -!)  +  (!  -  a)Ui(j) 

where  0<a<l . 

u-(j)  was  the  filtered  value  of  the  flow  u..  at  t,.  u^j)  was  the 

measurement  of  the  flow  u^  at  t. .  a=0.3  was  a  value  commonly  used. 

It  was  assumed  that  the  temperatures  were  measured  with  a  high  degree 

of  accuracy. 

III. 2. 2  Experimental  Work  for  Determining  6-|a 

The  experimental  runs  conducted  for  determining  6-|  and  92  were 
made  in  two  sets,  A  and  B.  The  runs  in  set  A  were  made  when  the  first 
effect  was  being  heated  and  9-.  was  estimated. 

Effects  1  and  2  were  initially  filled  to  their  steady  state 
hold-ups.  The  product  flow,  inter-effect  flow  and  feed  to  the  second 
effect  were  cut  off.  The  recirculation  flow,  W,-,,  was  held  between 
10  and  15  gpm  (0.6  to  0.9  kg/s).  Heating  steam  was  then  started  to 
the  first  effect.  Recordings  of  the  inter-effect  temperature  T-j2, 
steam  temperature,  T  ,  and  recirculation  flow,  W-i-j,  were  made  along 
with  the  other  variables  but  it  was  only  these  variables  that  figured 
in  the  ensuing  calculations.  Tables  3.1  to  3.5  contain  the  data  from 
Runs  Al  to  A5. 
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III. 2. 3  Calculations  and  Results  for  8, 

i  a 

A  computer  program  was  written  to  estimate  9,  from  the  data 
collected  in  Runs  Al  to  A5.  The  least-squares  program  used  was  sub- 
routine RMINSQ  (Westerberg,  1969).  This  program  was  based  on  a  program 
coded  by  M.  J.  D.  Powell  and  described  in  Powell  (1964).  This  routine 
has  the  capability  of  performing  a  least-squares  search  over  several 
functions  in  several  variables.  The  search  routine  does  not  require 
evaluation  of  derivatives. 

The  equation  describing  the  enthalpy  rise  in  the  first  effect 
is  equation  (3.3) . 

Runs  Al  to  A5  were  conducted  when  the  first  effect  was  not 
boiling  and  with  constant  hold-ups  in  the  first  and  second  effects. 
Thus,  the  liquid  entering  the  first  effect  was  only  the  recirculated 
liquid,  W,p  =  W-, ,  and  the  vaporization  was  zero,  V?,  =  0.  The  obser- 
vation T, 2  =  T-,  since  all  the  liquid  entering  the  first  effect  was 
recirculated.  For  every   value  of  9-,  which  subroutine  RMINSQ  searched 
over,  equation  (3.3)  was  integrated  from  the  initial  to  the  final 


ar"JrCMi2<hi2-hi>  +  v2i<hi  -  hP  +  <M  <3-3) 


2 
time  for  each  run.  Ten  functions  of  the  form  (T,   ,  -  J1      .     . ) 

v  1  ,calc    1 , observed' 

over  the  time  span  were  minimized  by  RMINSQ.  The  functions  toward  the 
end  of  the  time  interval  were  weighted  one  hundred  times  more  than  those 
at  the  start.  This  was  because  the  final  time  at  which  the  first  effect 
liquid  started  to  boil  was  more  important  from  the  point  of  view  of 
possible  switching  times  in  the  control  variables  at  this  point  in  time. 
The  results  of  the  minimization  are  tabulated  in  Tables  3.6  to  3.10  and 
these  values  are  plotted  in  Figures  3.2  to  3.6.  It  can  be  seen  that  the 
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calculated  values  are  indeed  least-squares  values  (with  the  weighting 

indicated)  which  converge  on  the  observed  values  toward  the  end  times. 

On  the  basis  of  these  results  a  value  for  6,  was  taken  as 

i  a 

0.1395. 

1 1 1. 2. 4  Experimental  Work  for  Determining  6-,,  and  6~ 

The  runs  in  Set  B  were  made  with  the  first  effect  boiling  and 
the  second  effect  initially  in  the  heating  and  later  in  the  boiling  stage. 
Three  parameters  were  actually  estimated  from  the  data  for  each  run. 
8,.  was  a  parameter  for  the  boiling  of  the  liquid  in  the  first  effect, 
9?  was  a  similar  parameter  for  the  heating  of  the  liquid  in  the  second 
effect  and  6?.  was  a  parameter  for  the  boiling  of  the  liquid  in  the 
second  effect.  Note  that  0-,  is,  in  effect,  the  revised  value  of  92a 
when  there  is  boiling  in  the  second  effect. 

The  experimental  procedure  consisted  of  bringing  the  hold-ups 
in  the  two  effects  to  their  steady  state  values.  The  controllers  in 
the  Transducing  and  Controlling  cabinet  (Chapter  II)  were  used  to  main- 
tain the  hold-ups  constant.  One  controller  was  used  on  the  first  effect 
and  another  on  the  second.  The  pneumatic  input  signal  to  the  first 
controller  was  the  output  from  the  DP  cell  which  measured  the  height 
(hold-up)  in  the  first  effect.  The  pneumatic  output  from  this  controller 
was  directed  to  valve  CV2  (Figure  2.2).  Thus  analog  control  of  the  first 
effect  hold-up  was  achieved  by  manipulating  flow  W-j2  which  is  the  liquid 
input  stream  to  the  first  effect  from  the  second  effect.  Likewise,  analog 
control  of  the  second  effect  hold-up  was  achieved  by  manipulating  flow 
WF  which  is  the  input  stream  to  the  second  effect. 

When  the  hold-ups  were  about  constant,  steam  was  let  in  to  the 
first  effect  for  heating  and  a  vacuum  of  around  ten  inches' mercury 
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(0.675  bar)  was  maintained  in  the  second  effect.  Data  were  recorded  when 
the  first  effect  temperature  was  near  boiling.  The  results  for  three 
runs  Bl,  B2  and  B3  are  shown  in  Tables  3.11,  3.12  and  3.13  respectively. 
III. 2. 5  Calculations  and  Results  for  elb  and  62 

A  similar  least-squares  search  using  subroutine  RMINSQ  was  used 
to  estimate  elb>  Qz&,   and  e2b  from  the  data  obtained  in  Runs  Bl  to  B3. 
The  function  evaluations  for  RMINSQ  entailed  integration  of  all  four  of 
the  differential  equations  (3.1)  to  (3.4)  as  neither  the  hold-ups  nor 
the  temperatures  were  held  constant.  Whenever  the  second  effect  tem- 
perature corresponded  to  the  temperature  of  boiling  in  the  second  effect 
(which  was  found  from  the  pressure  observed)  parameter  92b  was  used 
instead  of  B0   .  The  calculated  temperature  of  the  first  effect  solution 
T,  and  that  of  the  second  effect  solution  T2  were  obtained  from  the 
integration  of  the  state  equations  (3.1)  and  (3.3)  respectively.  The 
criteria  for  minimization  were  the  functions 

f(U-!l  <Tl,calc-Tl  .observed^  <3-3) 

n=l 


10 


f^=   .1   (T2,calc  "  T2.  observed^  (3'4) 

i=l 

Ten  functions  of  each  type  were  evaluated  in  the  time  span  of  each  run 
resulting  in  a  total  of  20  functions  for  the  evaluation  of  9lb>  62a  and 

92b- 

The  correspondence  between  the  observed  and  calculated  values  is 

shown  in  Tables  3.14  to  3.16,  while  Figures  3.7  to  3.9  are  plots  of  these 

values.  The  minimization  took  an  average  of  three  minutes  CPU  time  on 

the  IBM  370  for  each  run.  This  was  mainly  due  to  the  three-dimensional 
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search  and  the  integration  of  the  four  differential  equations  for 

every  function  evaluation.  Because  of  this  only  three  runs  were 

analyzed.  The  results  for  runs  Bl  and  B2  v/ere  much  better  than  those 

for  run  B3.  On  this  basis  the  mean  values  for  6-,,  ,  0O  and  0OU  were 

I  b   2a     2b 

taken  to  be  0.0928,  0.1359  and  0.1763  respectively.  Physically, 
this  meant  that  the  model  predicted  an  inside  film  coefficient  which 
was  from  80  to  90  percent  higher  than  that  obtained  experimentally. 
This  could  be  either  due  to  the  assumptions  made  in  the  model  or  to 
the  empirical  nature  of  the  film  coefficient  correlations  (A. 16),  (A. 39) 
and  (A. 51). 


CHAPTER  IV 
MINIMUM  TIME  CONTROL  POLICY 

This  chapter  deals  with  the  development  of  a  minimum  start- 
up time  control  policy  for  the  double  effect  evaporator  using  the  model 
equations  of  Chapter  III.  The  problem  is  stated  in  Section  IV. 1  and 
this  involves  identifying  state  and  control  variables,  equality 
constraints  in  the  form  of  algebraic  equations,  state  and  control 
variable  inequality  constraints,  and  possible  start-up  scenarios. 
Section  IV. 2  contains  the  derivation  of  a  general  algorithm  useful 
for  solving  minimum  time  problems  similar  to  that  for  the  evaporator. 
The  actual  use  of  the  algorithm  for  solving  the  present  problem  is 
described  in  Section  IV. 3.  It  contains  the  results  of  model  simulations 
in  arriving  at  the  optimal  policy  for  three  problems,  all  of  which  are 
minimum  time  start-up  problems  under  various  conditions.  Experimental 
verification  of  one  of  the  minimum  time  policies  and  the  effectiveness 
of  the  model  is  presented  in  Section  IV. 4. 

Refer  to  the  foldout  nomenclature  list  of  the  more  important 
symbols  at  the  end  of  this  chapter  (page  148)  to  aid  in  interpreting 
statements  made  using  these  symbols  in  the  succeeding  sections. 

IV. 1  Statement  of  the  Problem  for  the  Evaporator 

IV. 1.1  State  and  Control  Variables 

The  state  variables,  which  are  necessary  to  describe  completely 
the  state  of  the  process  at  any  particular  time,  are  the  differential 
variables  in  the  model  differential  equations  (3.1)  to  (3.6)  of 
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Section  III. 1.1.  For  uniformity,  x.  will  be  used  for  a  state  variable 
and  X  will  be  the  state  vector  with  components  x. .  These  are  assigned 
as  follows: 

x,  -  First  effect  hold-up  (I-!-,) 

x?  -  First  effect  liquid  enthalpy  (h-j) 

x3  -  Second  effect  hold-up  (M2) 

x4  -  Second  effect  liquid  enthalpy  (h2) 

xc  -  First  effect  solute  concentration  (C, ) 
5 

xr   -  Second  effect  solute  concentration  (C?). 

0 

The  decision  variables  listed  in  Section  III. 1.4  have  to  be 
defined  so  that  the  model  is  complete.  The  control  or  manipulative 
variables  are  chosen  from  this  set  depending  upon  the  controllability 
of  the  process  and  the  physical  realizability  of  the  control.  For 
example,  the  second  effect  vacuum  pressure  is  not  capable  of  being 
manipulated  physically  on  this  system  and  so  it  is  not  chosen  to  be  a 
control  variable.  The  feed  temperature  and  concentration  are  not  used 
as  control  variables  in  this  problem  either.  The  remaining  decision 
variables,  comprising  four  flow  rates  and  the  steam  pressure  to  the 
first  effect,  can  be  easily  manipulated  physically  and  can  be  used  to 
force  the  process  in  any  desired  direction.  For  example,  the  feed  to 
the  second  effect,  Wp,  and  inter-effect  flow  rate,  W-j2,  can  control 
the  inventories  in  the  first  and  second  effects.  The  steam  temperature, 
T  ,  and  recirculation  rate,  W,-|,  have  an  effect  on  the  first  effect 
temperature  and  the  rate  of  increase  of  the  second  effect  temperature. 
The  product  flow  rate,  kL, ,  is  used  to  control  the  product  concentration 
These  control  variables  are  capable  of  keeping  the  process  at  steady 
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state  and  the  steady  state  values  for  these  variables  are  governed  by 
the  steady  state  solution  of  the  differential  equations  (3.1)  to 
(3.6). 

For  uniformity,  let  u-  denote  a  control  variable  and  let  U  be 
the  control  vector  with  components  u-.  The  assignment  is  as  follows: 

u-,  -  Feed  to  the  second  effect  (W,-) 

u2  -  Intereffect  flow  rate  (W-jp) 

u3  -  Recirculation  flow  rate  (W-,-,) 

u,  -  Temperature  of  steam  to  first  effect  (T  ) 

iir  -  Product  flow  rate  out  of  first  effect  O'q-,). 

Rewriting  the  state  equations  (3.1)  to  (3.6)  and  the  algebraic 
equations  (3.7)  to  (3.12)  in  terms  of  the  state  and  control  variable 
nomenclature  defined  above,  we  have 


*1  =  W12  "  V21  "  u3  "  U5  (4J) 

X2  =  x~  [W12(h12  "  x2}  +  V21(x2  "  hl>  +  Ql]          (4-2) 

x3  =  ul  "  u2  '  V02  ^4"3^ 

X4  =  T   [ul(hF  "  X4}  +  V02(x4  "  h2}  +  Q2]            (4>4) 

V  ^12^12  "x5>  +V21X5^  (4-5) 

V^l^F"^  W  (4-6) 
Connection  equations, 

lJ12  =  u2  +  u3  ^"^ 

W12h12  =  U3X2  +  U2X4  (4.8) 
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W12C12  =  u3x5  +  u2x6  (4.9) 

Heat  transfer  equations, 

Q1  =  Q-j  (x-j ,  x2,  u3,  u4,  h12,  W12,  V21)  (4.10) 

Q2  -  Q2(x2,  x3,  x4,  hF)  (4.11) 

V21  =  Q2/X1  (4.12) 

IV. 1.2  State  and  Control  Variable  Constraints 

It  is  evident  that  in  a  real  system  the  control  variables 
cannot  take  on  all  values  as  there  are  physical  limitations  on  the 
maximum  and  minimum  flow  rates  and  temperatures.  The  lower  limit  for 
all  the  flow  variables  is  zero.  The  lower  limit  for  the  steam 
temperature  is  212°F  as  steam  cannot  be  supplied  at  a  lower  pressure 
than  atmospheric  in  the  first  effect.  The  upper  limit  depends  upon 
the  pipe  size  and  the  valve  size  for  the  flow  rates  and  on  the  steam 
supply  pressure  for  the  steam  temperature.  Thus,  all  the  control 
variables  are  subject  to  lower  and  upper  bounds  of  the  form 

u-   •  <  u-  <  u.  (4.13) 

i  ,mm  -  i  —  l  ,max  v    ' 

In  a  like  manner  some  of  the  state  variables  are  constrained. 
At  steady  state  all  the  state  variables  should  be  greater  than  or  equal 
to  their  steady  state  (desired)  values  (x.)--the  steady  state  hold-ups 
are  the  desired  operating  hold-ups,  the  steady  state  temperature  in 
the  first  effect  should  be  at  least  the  boiling  temperature  of  water  at 
1  atmosphere,  the  steady  state  temperature  in  the  second  effect  should 
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be  the  boiling  temperature  of  water  at  the  pressure  in  the  second 
effect,  and  the  steady  state  concentration  in  the  first  effect  should 
be  equal  to  the  desired  concentration.  The  upper  bounds  are  less 
clearly  defined;  for  example,  the  liquid  level  for  the  second  effect 
(hold-up)  should  not  exceed  the  overflow  limit.  The  upper  bounds 
on  the  temperatures  are  dictated  by  the  design  specifications  and  by 
characteristics  of  the  solution  being  concentrated.  In  general,  the 
state  constraints  are  given  by 

x1<x1  <x.>max  (4.14) 

IV. 1.3  Control  Scenarios 

In  the  start-up  of  the  evaporator,  it  is  useful  to  visualize 
the  change  with  time  of  the  state  variables  for  certain  values  of  the 
control  variables.  It  is  evident  that  the  control  variables  determine 
the  order  in  which  the  state  variables  reach  their  steady  state  or 
desired  values.  Intuitively,  the  optimal  policy  will  endeavor  to 
force  each  state  variable  directly  to  its  final  steady  state  value 
and  then  maintain  it  while  bringing  the  others  to  steady  state.  The 
order  in  which  the  state  variables  reach  their  steady  state  values 
completes  the  scenario.  Theoretically,  there  should  be  a  total  of  n 
factorial  scenarios  for  a  system  with  n  state  variables.  But  most  of 
these  are  not  possible  as  certain  state  variables  can  reach  steady 
state  only  if  certain  others  have.  For  example,  in  the  case  of  the 
evaporator,  the  first  effect  liquid  has  to  boil  before  the  second 
effect  liquid  does.  The  equations  describing  a  stage  in  a  scenario 
are  different  from  those  describing  another  stage.  For  each  stage  they 


68 


are  simplifications  of  the  general  equations.  A  typical  scenario 
for  the  start-up  of  the  evaporator  (and  the  resulting  simplifications 
in  the  general  equations  (4.1)  to  (4.12))  is  shown  below. 
Stage  A:  tn  <_  time  t  <  t-,  .  Filling  and  heating  first  effect. 
Control  variables: 

Feed  to  second  effect  at  maximum  flow  rate,  u-,  =  u-. 
All  feed  delivered  directly  to  first  effect,  u9  =  u-. 
Temperature  of  steam  for  first  effect  at  maximum  value,  u»  =  u* 
No  recirculation  or  product  flow  possible,  u^  =  u5  =  0. 
Resulting  state  equations: 

First  effect  hold-up  increasing;  x-,  =  u2,  x-iUq)  =  0 

Heating  of  first  effect;  x2  =  —  [u2(x4  -  x2)  +  Qi]»  x2(tQ)  =  lv 

No  increase  in  second  effect  hold-up;  x~  =  0,  x-  =  Xo(0  =  0 

No  heating  of  second  effect;  x»  =  0,  x.   =  x»(t«)  =  hp 

No  concentration  occurring  in  either  effect;  kr  =  kr  =   0; 

X5  =  X6  =  X5(t0}  =  X6(t0}  =  CF- 
Time  t-, ,  signifying  the  end  of  Stage  A,  is  determined  when  the  first 

effect  is  filled,  i.e.,  when  x,  (t-,)  =  x, . 

Stage  B:  t-,  <_  time  t  <  t?.  First  effect  solution  is  being  heated 

and  second  effect  is  being  filled. 

Control  variables: 

Feed  flow  to  second  effect  at  maximum,  u^   =  un 

1    1 ,max 

Flow  to  first  effect  stopped,  iu  =  0 

Recirculation  flow  set  at  maximum,  u0  =  u0  „„ 

3    3, max 

Temperature  of  steam  to  first  effect  set  at  maximum  value, 
u4  =  u4,max 
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No  product  withdrawn,  Ur  =  0. 
Resulting  state  equations: 

No  change  in  first  effect  hold-up;  x-,  =  0,  x-,  =  x-.  (t-. )  =  x-. 
Heating  of  first  effect;  x?  =  Q-i/x-, 
Second  effect  hold-up  increasing;  x3  =  u, ,  x3(t, )  =  0 
No  heating  of  second  effect;  x.  =  0,  x.  =  x,(t-,)  =  lv 
No  concentration  occurring  in  either  effect;  x5  =  xg  =  0, 

X5  =  X6  =  X5(V  =  X6(V  =  CF- 
Time  t?,  signifying  the  end  of  Stage  B,  is  determined  when  the  second 

effect  is  filled,  i.e.,  when  x3(t2)  =  x3< 

Stage  C:  t2  £  time  t  <  t3<  First  effect  solution  is  heated  to  boiling, 

Control  variables: 

Feed  flow  stopped,  u-,  =  0 

Flow  to  first  effect  stopped,  u2  =  0 

Recirculation  flow  set  at  maximum,  u0  =  u0  „  „ 

3    3, max 

Temperature  of  steam  to  first  effect  set  at  maximum  value, 

u4  =  u4,max 
No  product  withdrawn,  Ur  =  0. 

Resulting  state  equations: 

No  change  in  first  effect  hold-up;  x,  =  0,  x,  =  x-,(t2)  =  x. 

First  effect  is  being  heated;  x?  =  Q-./X-. 

No  change  in  second  effect  hold-up;  x3  =  0,  x3  =  x3(t2)  =  x3 

No  heating  of  second  effect;  x»  =  0,  x«  =  x*(t«)  =  hp 

No  concentration  occurring  in  either  effect;  x5  =  Xg  =  0; 

x5  =  x6  =  X5(t2}  =  x6(t2>  =  CF- 
Time  t~s   signifying  the  end  of  Stage  C,  is  determined  when  the  first 

effect  starts  to  boil,  i.e.,  when  x2(t3)  =  Xp« 
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Stage  D:  t-  <  time  t  <  t4-  First  effect  solution  is  boiling  (and 
becoming  concentrated).  Second  effect  solution  is  being  heated. 
Control  variables: 

Feed  to  first  effect  is  equal  to  vaporization  from  first 
effect  to  maintain  constant  hold-up,  u2  =  V21 

Feed  to  second  effect  set  to  maintain  constant  hold-up  in 
second  effect,  u-,  =  u^ 

Recirculation  at  maximum,  u^  =  u^  max 

Temperature  of  steam  to  first  effect  at  maximum,  u4  =  Uqtmx 

No  product  withdrawn,  Ur  =  0. 
Resulting  state  equations: 

No  change  in  first  and  second  effect  hold-ups;  x-j  =  X3  =  0; 

xl  =  x1^3^  =  xl '  x2  =  x2^3^  =  X2 

First  effect  is  being  heated;  x«  =  —  [u2x4  +  Q]  -  V21h^], 

xl 
X2  3  ~~  x2 

Second  effect  is  being  heated;  x4  =  7-  [u-,(hF  -  x^)  +  Q2] ; 

x3 
x4(t3)  ■  hp 

First  effect  solution  is  being  concentrated;  x5  =  —  (u2Cp); 


xl 
x5(t3)  =  CF 

No  concentration  of  second  effect  solution;  xg  =  0;  xg  =  Xg(t3)  =  Cp, 
Time  t*   signifying  the  end  of  Stage  D  is  determined  when  the  second 
effect  starts  to  boil,  i.e.,  when  x4(t4)  =  x4> 

Stage  E:  t4  <_  time  t  <  tf.  Solution  in  both  effects  is  boiling  and 
being  concentrated. 
Control  variables: 

Feed  to  first  effect  set  equal  to  vaporization  in  first 
effect  to  maintain  constant  hold-up,  u2  =  V21 
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Feed  to  second  effect  set  to  maintain  constant  hold-up  in 

second  effect,  u,  =  u?  +  Vfi? 

Recirculation  at  maximum,  u0  =  u0 

3    3, max 

Temperature  of  steam  to  first  effect  at  maximum,  u»  =  u. 
No  product  withdrawn,  Ur  =  0. 
Resulting  state  equations: 

No  change  in  first  or  second  effect  hold-ups;  x-,  =  x3  =  0; 

X-,  -  X-j^t»J  -  X-j  ,  Xo    Xo\t^j  -  x. 

First  effect  is  being  heated;  x?  =  —  (u?x.  +  0,  -  V9,h,) 

xl. 
No  change  in  second  effect  enthalpy;  x«  =  0,  x»  =  x«(t.)  =  x. 

First  effect  solution  is  being  concentrated;  x5  =  —  (u-iXg) 

xl 
Second  effect  solution  is  being  concentrated; 

x6  =  —  (u^p  -  u2x6),  x6(t4)  =  Cp. 

X3 
Time  tf,  signifying  the  end  of  Stage  E,  and  consequently  the  final 

time,  is  determined  when  the  first  effect  concentration  reaches  a 

desired  value,  i.e.,  when  Xr(tf)  =  Xr. 

After  the  final  time,  tf  found  above,  the  control  variables 

u-, ,  u2  and  Ur  are  obtained  from  the  steady  state  solutions  of  the 

differential  equations  (4.1),  (4.3)  and  (4.5).  The  feed  to  the  first 

effect  is  such  that  a  constant  hold-up  is  maintained  in  the  first 

effect,  i.e.,  \x?   =  Ur  +  V?, .  The  feed  to  the  second  effect  is  such  that 

a  constant  hold-up  is  maintained  in  the  second  effect,  i.e.,  u-,  =  u^  +  V^' 

The  product  rate  is  such  that  a  constant  product  concentration  is 

x6 
obtained,  i.e.,  Ur  =  u0   —  .  Simplifying  the  above  three  relationships 
5    d   x5 

we  can  obtain  a  sequence  for  determining  the  control  variables  u, , 
u?  and  Ur  at  steady  state  as  follows: 
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u2  =  V21/(l  +^]  ;  U]  =u2  +  VQ2;  u5  =  u2^ 

Note  that  certain  point  constraints  of  the  form  x^t-j)  =  xi 
separate  the  various  stages  of  the  scenario.  The  control  variables, 
state  and  algebraic  equations  change  when  these  points  in  time  are 
encountered. 

It  is  assumed  that  certain  state  variables  such  as  the  first 
and  second  effect  hold-ups  will  be  maintained  at  their  steady  state 
values  once  these  are  attained.  That  is,  the  point  constraint 
x.(t.)  =  x.  is  followed  by  the  equality  constraint  x^  =  0.  This 
latter  equality  is  maintained  by  calculating  the  required  value  of  an 
appropriate  control  variable.  This  is  similar  to  the  concept  of  a 
first  order  state  variable  inequality  constraint  (Bryson  et  al . ,  1963), 
(Bryson  and  Ho,  1969). 

It  is  clear  that  another  scenario  will  give  rise  to  a  different 
ordering  of  the  point  constraints,  x^t.)  =  x.s  depending  upon  the 
order  in  which  the  state  variables  reach  their  steady  states.  The 
optimal  scenario  is  the  one  which  will  result  in  a  minimum  final  time. 

This  scenario  approach  simplifies  the  mathematical  problem 
to  a  great  extent  as  a  minimum  number  of  state  equations  have  to  be 
integrated.  The  equations  simplify  considerably  leading  to  simpler 
adjoint  equations  and  Hamiltonian  minimizations. 

IV. 1.4  Summary  of  the  Problem  Statement 

The  problem  is  to  minimize  the  final  time,  that  is 

Min  J  =  tf 
U 
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subject  to  the  state  equations  (4.1)  to  (4.6),  the  connection  equations 
(4.7)  to  (4.9),  and  the  heat  transfer  equations  (A.l)  to  (A. 52)  where 
the  control  variables  are  constrained, 


U  .  <  U  <  U 
mm  —  —  max 


and  point  constraints  pertaining  to  the  particular  scenario  are  to  be 
satisfied 

x.j(t,0  =  *j  ;   i  =  l,n 

J  =  1  ,n 

IV. 2  A  Minimum  Time  Algorithm 

IV. 2.1  General  Problem 

The  objective  is  to  minimize  the  final  time  by  selection  of 
the  controls  U 


Min  J  =  t, 
U      1 


subject  to: 

the  state  equations 


(4.15) 


X  =  f(X,  U,  Z)  ;  x(0)  given 
the  algebraic  equations 

g(X,  U,  Z)  =  0 
the  control  constraints 


u  .  <  u  <  u 

mm  -   —  max 


or 


hli(ui}  -   ui  "  ui,max^ 


h2i(ui}  =  ui,min  -  ui  ^° 


h(U)  <  0 


[4.16) 


[4.17) 


(4.18) 


(4.19) 
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and  the  point  constraints 


x0(t0)   xQ 


xkl^kl^   xkl 

xk2(tk2)  =  xk2  (4.20) 

t„  =  fixed 

IV. 2. 2  Lagrange  Formulation  and  Necessary  Conditions 

Assume  that  the  point  constraints  are  met  at  times  t«,  t, ,,..., 
tj;.  Also  assume  the  times  t.  correspond  to  all  those  times  when  one 
(or  more)  of  the  controls  or  states  reaches  or  leaves  a  constraint. 
Next,  form  the  following  index  sets, 

I0  -  {0,  1,  2,  kl  ,  k2 ,  f} 

l}   =  {1,  2 kl,  k2 ,  f} 

I2  =  {0,  1,  ,  kl,  k2,  f-1} 

Within  each  set,  order  the  indices  so  the  times  t.,  i  £  L  or  I,  or  I« 
are  in  increasing  order.  Also  form  the  index  set 

K  =  {kl,  k2 } 

Introducing  multipliers  a,  6  and  A  we  can  write  the  Lagrangian 

for  the  problem 

L(X,  U,  t.,  i  £  IQ,  a,  6,  X)  =   tf  +  af(xf  -  xf(tf))  +  I     ak(xk  -  xk(tR)) 
t"  ^e^ 

+  I  (>T(f  -  x)  +  BTh)dt  (4.21) 

ieli  + 


1  'i-i 
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The  algebraic  equations  g(X,  U,  Z)  =  0  are  not  included. 
These  will  always  be  satisfied  by  solving  for  the  dependent  variables 
Z  and  substituting  the  resulting  values  into  the  state  equations. 
Also  the  equations  at  t~,  X(tfi)  =  X~,  will  always  be  satisfied  and  are 
not  included.  Note  that  the  fourth  term  accounts  for  the  changes  in 
the  state  equations,  algebraic  equations,  and  control  variable 
constraints  along  the  trajectory.  This  term  has  been  written  to  allow 
for  possible  discontinuities  in  the  Lagrange  multipliers,  A,  at  the 
entry  and  exit  corners  of  constraints. 

We  define  the  Hamiltonian  H  =  A  f  and  rewrite  equation  (4.21) 

L(X,  U,  t.,   i  e  IQ,  a,  3,  A)  =  tf  +  af(xf  -  xf(tf)) 

f*1      T    T 

+  I    a.(x.  -  x.(t.))  +  I  (H  -  A'X  +  3'h)dt     (4.22) 

keK  K  K    K  K     iel1j  + 

Ti-1 

On  an  extremal  solution,  the  Lagrangian  L  must  be  stationary 
with  respect  to  small,  arbitrary  perturbations  of  the  times  t.(i  e  I,); 
the  multipliers  a,  6,  and  A;  and  the  states  X  and  controls  U.  To 
derive  the  necessary  conditions  for  a  stationary  point  of  the  Lagrangian, 
we  take  first  order  variations  of  L  with  respect  to  t. (i  e  I-i),  a,  3> 
A,  X  and  U, 

6L  =  (1  +  H  -  ATX  +  3Th  -  afxf)t  6tf  -  af6xf(tf)  +  (xf  -  xf(tf))6af 

+  kU~  ak5Xk(tk}  +   ("k  "  MVK  "   (Vk}tk6tk^ 

+     I       [(H  -  ATX  +  3Th).-6t.    -   (H  -  ATX  +  3Th).+     6t.    ,] 
iel2  ri     ^  ri-l     1_1 
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s- 


1    t+ 
l-l 


+     l  .      f|H  _x]'    6A+    [3H     +gT3h     _AT6x]5X+    [  3H 

ill,   J        I    i  3X        J  i  3XT  8XT  J  I  3UT 


+  3T  ^V  IfiU  +  hT66 
3U1 


(4.23) 


Integrating  the  term 


A  SXdt  by  parts, 


'1-1 


i 


1 


AT6xdt  =  (AT5x)  +   -  (AT6x)  _  + 

+  ti-1       *1   t+ 

ri-l  ri-l 


AT6xdt 


[4.24' 


Substituting  (4.24)  into  (4.23)  and  collecting  coefficients  of 
6a,  66,  6A,  6X,  6U,  6x,  and  6t,  separately,  we  have 

TV,  .  J, 


6L  =  (1  +  H  -  AX  +  B  h  -  afxf)t  6tf 
-  af6xf(tf)  -  AT(tf)6X(tf) 
+  AT(t0)6X(tQ) 
+  (xf  -  xf(tf))6af 


+  I       (AT)  +6X(t.)  -  (AT)  6X(t.) 
ieL 

i^K 


lei2   t: 


ti 


(a) 
(b) 
(c) 
(d) 
(e) 


Tc  .  „T, 


+  I       [(H  -  ATX  +  6Th)  6t.  -  (H  -  X'X  +  e'h)  ,6t. 

id2  t:  t| 


f] 


+  I  [(\T)  +6X(tk)  -  (AT)  _6X(tk)  +  ak6xk(tk)J    (g] 


keK 
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+  I     (xk  "  xk(tk))6ak  (h) 

keK 

+  I     [(H  -  ATX  +  6Th)  6t,  -  (H  -  ATX  +  BTh)  ,6t, 
keK  t"  tk 

+  (akXR)t  6tk]  (o) 

MtT    (ix-*)Tsxdt  (p) 

lei-,  . 

ri-l 

t- 

+  I       [AT  +  ^V+  3T^V]6xdt  (q) 

1eI,J  +        3X1      9X1 
1  ti-l 

t- 
+  I  [  3"  +  gT  3h  ]fiUdt  (r) 

iel-,  J+     3U1      3U1 
ti-l 

^       T 
+  I  h'63dt  (s) 

ieIi  ;+ 

i~1  (4.25) 

The  Kuhn-Tucker  conditions  arise  from  the  Kuhn-Tucker  multipliers  3  and 
the  inequality  constraints  h, 

eohj  ■  °  (t) 

B>0  (u) 

For  stationarity  of  the  Lagrangian  all  the  coefficients  of  6X, 
6U,  6x,  ,  6a,  63,  6X  and  6t-  occurring  in  equation  (4.25)  must  be  zero. 
By  equating  these  coefficients  to  zero  we  arrive  at  the  following 
necessary  conditions  (NC). 
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From  term 

From  term 

From  term 

From  term 

From  term 

From  term 
From  term 
From  term 

From  term 
From  term 

From  term 

From  term 
From  term 
From  term 


q) ;  X j   -  3  — f  ;  tj  ,  <  t  <  t. 

dr  9X1    1_l       ] 

r);  i"+  6T^V=  0;  tt   <  t  <  tT 
3U1      9U1       1_l        1 

p);  X  =  |£=  f(X,  U);  t|_1  <  t  <  tT 


t);  3-h-  =  0;  tj_1  <  t  <  tT 


u);  6  >  0;  t.  ■,  <  t  <_  t". 


h);  xk  -  xk(tk)  =  0;  ke  K 


d);  xf  -  xf(tf)  =  0 


b);  Af(tf)  =  -  af(tf) 

Xj(tf)  =  0;  j=l,...,n;  j7f 


c);  ^i(^n)  =  unknown  since  6x(t0)  =  0 

e);  Aj (tt)  =  A^tT);  j=l n 

i  e  1 2 »  H& 


g);  X.(tk)  =   A^t");  j=l,...,nj  j/k;  k£K 
Ak(tk)  =  Ak(tk)  -  ak;  keK 


a),  (NC4),  (NC5)  and  (NC8);  H(tf)  =  -1 

f),  (NC4)  and  (NC11);  H(tT)  =  H(t|);  ieI2>  i£K 

o),  (NC4)  and  (NC12);  H(t~)  =  H(t£);  keK 


(NCI) 

(NC2) 

(NC3) 

(NC4) 

(MC5) 

(NC6) 
(NC7) 

(NC8) 
(NC9) 

(NC10) 

(ricn) 

(NC12) 
(NC13) 
(NC14) 
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IV. 2. 3  Comments  on  the  Necessary  Conditions 

The  state  equations  and  point  constraints  are  the  necessary 
conditions  NC3,  NC6  and  NC7.  The  Kuhn-Tucker  conditions  NC4  and 
NC5  indicate  that,  when  a  control  constraint  h  is  encountered,  that 
constraint  is  held  by  choosing  U  such  that  h(X,  U)  =  0  provided  that 
the  Kuhn-Tucker  multiplier  B  is  non-negative. 

The  problem  is  set  up  such  that  the  point  constraints  arise 
when  there  is  a  change  in  the  form  of  the  state  equations  and, 


consequently,  a  change  in  the  adjoint  equations  and  Hamiltonian. 
For  example,  before  tim 
the  state  equations,  or 


For  example,  before  time  t,  (where  x-^tj)  =  x-| )  let  f^   represent 


x  =  f(])(x,  u)       t  <  t1 

After  time  t,  let  f^  '   represent  the  state  equations,  that  is 

X  =  f(2)(X,  U)         t  >  t1 
The  Hamiltonians  are 

hM-xV1'  and  H<2'=ATf<2> 

The  functional  forms  f^  and  V2'   reflect  the  change  in  the  state 
equations.  Condition  NCI 4  shows  that  the  Hamiltonian  is  continuous 
across  the  times  t,  where  point  constraints  are  encountered.  For 
example,  at  t  =  t-j ,  x-^t-j)  =  x-j  and  \v   '  =  \v     . 

Condition  NC11  states  that  the  Lagrange  multiplier  (or  adjoint 
variable),  corresponding  to  the  particular  state  on  which  there  is  a 
point  constraint  at  t.,  has  a  discontinuity  at  t^ ,  whereas  the 
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multipliers  corresponding  to  the  other  states  are  continuous  at  ty 
At  time  t.,  for  example, 

X-,(t|)  =  X-^tf)  -  a-, 

X.(tj)  =  X..(t~)  ;  j=2 n 

The  value  of  the  "jump,"  ^  ,  in  the  multiplier  A-,,  is  readily  determined 
from  the  condition  H^  '  =  \r    K 

At  times  t.,  i  e   IQ,  when  control  constraints  have  to  be 
satisfied,  all  the  multipliers  are  continuous  as  the  trajectories 
enter  and  leave  the  constraints— NC10.  The  Hamiltonian  is  also 
continuous  at  such  points  in  time  as  shown  by  condition  NCI  3. 

flecessary  conditions  NC8  and  NC12  provide  values  for  the  adjoint 
variables  and  Hamiltonian  at  the  final  time.  The  adjoint  variables 
corresponding  to  the  states  that  are  unconstrained  at  the  final  time 
have  a  value  zero.  If  only  one  state  is  constrained  at  the  final 
time,  the  adjoint  variable  corresponding  to  this  state  can  be 
determined  from  the  final  condition  on  the  Hamiltonian 

H(tf)  =  -1 

However,  if  more  than  one  state  variable  is  constrained  at  the  final 
time,  the  values  for  all  but  one  adjoint  variable  have  to  be  guessed 
and  the  remaining  one  adjoint  variable  can  then  be  determined  from  the 
final  condition  on  the  Hamiltonian.  Let  us  suppose  that  the  adjoint 
variable  corresponding  to  the  stopping  condition  x^(tf)  =  xf  is 
determined  from  the  final  value  of  the  Hamiltonian  and  that  another 
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variable  A.  has  been  guessed.  On  subsequent  iterations  this  value  of 
X.   is  updated  by  noting  that  the  gradient  of  L  with  respect  to  X.(tf) 
is  -  (x.  -  x.(tp)).  The  Saddle  Point  Theorem  requires  us  to  maximize 
L  with  respect  to  X-(tf)  which  is  equivalent  to  driving  the  gradient  to 
zero.  So  the  value  of  X.(t-)  on  following  iterations  should  be  in  such 
a  direction  that  at  the  final  time  the  deviation  of  x. (t-)  from  the 
desired  value  x. ,  i.e.,  (x.  -  x. (t-))  is  driven  to  zero. 

This  same  iterative  technique  is  used  if  more  than  one  point 
constraint  is  encountered  at  an  intermediate  point  in  time  t.  .  One 
of  the  adjoint  variables  at  time  tj"  can  be  determined  by  using  the 
continuity  of  the  Hamiltonian  at  time  t, .  The  other  variables  will 
have  to  be  guessed  initially  and  then  updated  using  the  gradients 
available  on  subsequent  iterations. 

Necessary  condition  NC2,  along  with  the  Kuhn-Tucker  conditions 
NC4  and  NC5  enables  us  to  replace  condition  NC2  by  another  condition 
made  possible  by  invoking  the  strong  minimum  principle.  Denn  (1969) 
has  lucidly  shown  how  the  Hamiltonian  takes  on  its  minimum  value  for 
the  optimal  decision  function  U(t),  both  on  and  off  the  inequality 
constraints.  Utilizing  this  result,  we  can  replace  the  stationary 
condition  NC2  by 

Min  H(X,  U,  X) 
U 

subject  to  U  .  <  U  <  U 
J      mm  -  -  max 

IV. 2. 4  Minimum  Time  Algorithm 

An  algorithm  for  arriving  at  the  minimum  time  policy  was 
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developed  based  on  the  necessary  conditions  arrived  at  in  Section 
IV. 2. 2.  It  can  be  classified  as  a  "Minimum  H"  algorithm  as  it  deals 
with  direct  Hamiltonian  minimizations  along  the  trajectory.  Proposing 
a  control  scenario  to  start  with  is  essential  to  simplify  the  problem 
and  to  use  the  algorithm  effectively. 

Let  the  various  elements  of  the  state  vector  be  divided  into 
3  groups  as  follows: 

Group  A:  States  which  remain  unconstrained  through  start-up. 

Group  B:  States  which  meet  their  steady  state  values  during 
start-up  and  which  define  the  point  constraints. 

Group  C:  States  which  meet  their  steady  state  values  at  the 
final  time. 

The  state  equations  are  modified  along  the  trajectory  when 
the  point  constraint  conditions  are  met  by  the  variables  in  Group  B. 
Also,  all  algebraic  equations  are  satisfied  throughout  the  trajectory, 
and  this  is  implied  in  the  state  equations.  When  a  variable  belonging 
to  Group  B  arrives  at  its  steady  state  value,  it  is  assumed  that  for 
subsequent  time  the  state  equation  is  replaced  by  the  corresponding 
algebraic  equation  x  =  0.  Thus  the  Hamiltonian  H  =  X   f  is  different 
for  points  on  the  trajectory  depending  upon  which  of  the  state 
equations  are  active;  so  also  are  the  adjoint  equations 

A  -  -  — y  -  g  — j     ,     t.  ,  <_  t  <_  t. 
3X1      3X'      1_l       1 

The  algorithm  proceeds  as  follows: 

1)  Guess  a  nominal  control  policy  U  which  will  cause  a  stopping 

condition,  x.  =  xf,  to  be  satisfied  at  the  final  time. 
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2)  Integrate  the  state  equations  forward  until  the  stopping 
condition  is  satisfied.  This  implicitly  determines  the  final  time. 
The  proper  set  of  state  equations  should  be  integrated  depending  upon 
the  point  constraints,  the  control  constraints,  and  the  algebraic 
equations  which  have  to  be  satisfied  along  the  trajectory. 

3)  At  the  final  time  determined  in  step  (2),  let  all  the 
multipliers  of  the  variables  in  Group  A  be  zero.  Guess  the  multipliers 
of  all  the  variables  in  Group  C  with  the  exception  of  that  corresponding 
to  the  stopping  condition.  Determine  this  latter  multiplier  from  the 
final  condition  on  the  Hamiltonian, 

Vf =  -  ]  -  .L  Vi 
w 

4)  With  the  values  of  the  adjoint  variables  at  the  final  time 
determined  in  step  (3),  integrate  the  adjoint  equations  in  the  reverse 
direction.  At  times  t,  ,  when  point  constraints  were  met  on  the  forward 
integration,  determine  the  values  A(tj")  by  utilizing  the  continuity 
and  jump  conditions  on  the  adjoint  variables  and  Hamiltonian.  For 
example,  if  only  one  point  constraint  of  the  form  x, (t. )  =  x.  is  met 

at  time  t.  ,  then  ^(tj")  can  be  determined  from  the  continuity  of  the 
Hamiltonian 

J7k  j/k 

and  A.(tj")  =  A.(t.  );  jen,  j^k.  If  more  than  one  point  constraint  is 
met  at  time  t.  ,  then  the  values  of  all  but  one  of  these  multipliers 
should  be  guessed  at  time  tj~  and  the  last  one  determined  from  the 
continuity  of  the  Hamiltonian. 
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5)  Simultaneously  on  the  reverse  iteration,  minimize  the 

T  * 

Hamiltonian  H  =  A  f  at  each  point  to  determine  the  optimal  U 

Min  H(X,  U,  A) 

* 
U 

6)  On  reaching  time  tQ,  update  the  control  policy  used  on 
the  forward  integration  with  that  found  on  the  reverse  integration  in 
step  (5) 

Ui+1  =  U1'  +  c(t)(U*  -  U1') 

c(t)  is  chosen  to  limit  the  change  in  U  if  too  large  a  change  is 
indicated. 

7)  Integrate  the  state  equations  forward  as  in  step  (2). 
At  the  final  time  determine  the  difference  in  the  states  for  the 
variables  in  Group  C  from  their  desired  values  and  update  the  guess  on 
the  multipliers  A.  such  that  the  gradient  of  L  with  respect  to 
A.(tf)  is  driven  to  zero,  i.e.,  to  drive  (x.  -  x.(t.))  to  zero.  As 
before,  determine  the  A  corresponding  to  the  stopping  condition  from 
the  final  value  of  the  Hamiltonian. 

8)  Integrate  the  adjoint  equations  in  the  reverse  direction 
as  in  step  (4).  Update  the  guess  on  A.(tJ~);  jen,  j7k  in  a  similar 
manner  as  in  step  (7).  Determine  A,  as  before  from  the  continuity 

of  the  Hamiltonian. 

9)  Repeat  steps  (5)-(8)  until 

(1)  6J  <  e-j  -  no  significant  improvement  in  the  final  time 

(2)  IIU1       -  U1  W<  £p  -  no  significant  change  in  the  control   policy 

and  (3)       II6A.II     <  n  at  tf 

and  (4)        II6A.H     <  n  at  t" 

J  K 
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The  optimal  policy  is  chosen  as  the  one  which  satisfies  the 
above  conditions. 

IV .3  Solution  to  the  Evaporator  Problem 
IV. 3.1  Problem  1.  Constraint  on  the  Second  Effect  Hold-up 

The  problem  was  solved  for  the  scenario  described  in 
Section  IV. 1.3.  The  concentration  equations  (4.5)  and  (4.6)  were  not 
used  in  this  simulation.  The  Hamiltonians  and  adjoint  equations  for  the 
various  stages  of  the  scenario  are  given  below. 

The  general  Hamiltonian  for  the  problem  is: 

H  =  Xl(W12  -  V21  -  u3)  +  ^  [W12(h12  -  x2)  +  V21(x2  -  h{)  +  Q]] 

+  X3(Ul  -  u2  -  VQ2)  +  ji  [Ul(hp  -  x4)  +  VQ2(x4  -  hp  +  Q2] 

This  simplifies  for  the  various  stages  as  follows: 

Stage  A:  tn  <  t  <  t, 

A2 
Ha  =  A1U2  +  ~  ^2^4  "  x2^  +  *V  +  MU1  "  u2^ 


Stage  B:  t,  <  t  <  t?.  x, (t,)  =  x, 

Hb  =  ^l  +  A3(ul  -U2} 
xl 

H|2  =  u«  ;  u~  =  0  or  f,  =  0 
Stage  C:  t?  <  t  <  t,     x-,(t?)  =  x,, 

Xl 

u-,  =  u2  or  f~  =  0  and  f -,  =  0 
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Stage  D:     t~  <  t  <  tf  x?(t^)  >  x? 

X? 
Hd  =  ~[W12(h12  -  x2)  +  V21(x2  -  hp  +  Q^ 

xl 

+  ^  [Ul(hF  -  x4)  +  q2] 

x3 
W-J2  =  Vp-i  +  u2  or  f,  =  0;  u-,  =  u2  or  f_  =  0 

Time  t.  is  determined  when  the  second  effect  starts  to  boil,  i.e., 
when  x»(tf)  =  x». 

The  variables  Q-, ,  Q2,  V^-, ,  VQ2,  W-,,,,  h,,,,  h^,  lv  are  determined 
from  the  algebraic  equations  in  Appendix  A. 

The  adjoint  equations  are 

and  H  is  either  H  ,  H,  ,  H  or  H,. 
a   b   c    d 

9H 
The  partial  derivatives  -rr  are  evaluated  numerically  for  each  point 

in  time  using  the  appropriate  Hamiltonian  valid  at  that  time.  Note 

that  in  Stage  A,  X^  =  0;  this  is  because  the  liquid  enthalpy  in  the 

second  stage,  x.,  is  constrained  to  stay  at  x.  =  !v  by  xA  =  0.  In 

Stage  B,  X-j  =  0  since  the  first  effect  hold-up  X-.  =  X-,  =  constant. 

That  is,  the  objective  function  is  made  insensitive  to  changes  in  x-, 

by  holding  x1  =  x-j .  Likewise,  in  Stage  C,  X3  =  0  and  X1  =  0  because 

both  holdups  are  fixed. 

The  final  time  is  determined  when  the  second  effect  boils,  i.e., 

when  x4(tf)  =  x^.  At  tf,  X2  and  X4  are  active.  X2  =  0  since  the  final 

condition  on  the  liquid  enthalpy  in  Stage  1,  x0,  is  unspecified.  X.  is 
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determined  from  the  final  condition  on  the  Hamiltonian.  On  the 
reverse  integration  of  the  adjoint  variables  at  time  t«  the  second 
effect  holdup,  x3,  becomes  constrained.  The  mulitplier  Xo(tZ)  is 
found  from  continuity  of  the  Hamiltonian  H.  (t«)  =  H  (t2).  At  time 
t-, ,  the  first  effect  holdup,  x-, ,  becomes  constrained  and  A-,(t7)  is 
determined  from  the  continuity  of  the  Hamiltonian  at  time  t-, , 
H  (t7)  =  H,  (t-,).  Also,  on  the  reverse  integration,  the  appropriate 
Hamiltonian  for  each  stage  is  minimized  with  respect  to  the  control 
variables  U.  Note  that  in  Stage  A  minimization  of  H  is  with  respect 
to  second  effect  feed  u-, ,  intereffect  flow  u?,  and  first  effect  steam 
temperature  u*.  The  recirculation  rate  u,  and  the  product  flow  u5 
cannot  be  started  when  the  first  effect  is  still  filling.  In  Stage  3 
the  feed  flow  to  the  first  effect,  Up,  has  a  set  value  as  it  maintains 
the  first  effect  holdup  x-,  =  x-, .  So  minimization  of  H,  is  with  respect 
to  Up >  u3  and  u,.  Likewise,  in  Stage  C,  the  second  effect  feed  flow, 
u-, ,  is  fixed  to  maintain  the  second  effect  holdup  at  x,  =  x~  and 
minimization  of  H  is  with  respect  to  the  recirculation  flow  u,  and 
the  steam  temperature  u».  The  same  holds  for  Stage  D. 

Minimization  of  the  Hamiltonian  was  achieved  by  a  slightly 
modified  version  of  the  computer  program  VA04A  originally  coded  by 
M.  J.  D.  Powell.  It  is  based  on  a  conjugate  gradient  search  technique 
which  does  not  involve  partial  derivatives,  and  it  is  explained  in 
Fletcher  and  Powell  (1963).  This  routine  was  adapted  for  bounded 
variable  minimizations  by  writing  a  package  which  accounted  for  the 
control  variable  constraints.  When  a  control  constraint  was  encountered, 
a  perturbation  of  the  control  variable  into  the  feasible  region  decides 


whether  the  constraint  should  be  held  or  be  released.  This  is 

basically  a  numerical  evaluation  of  the  associated  Kuhn-Tucker  multiplier. 

The  forward  and  reverse  integrations  were  done  using  Hamming's 
predictor-corrector  method  HPCG  available  in  the  IBM  Scientific 
Subroutine  Package.  The  adjoint  equations  required  the  numerical  values 
of  partial  derivatives,  and  a  subroutine  PDERIV  was  written  to  do  this. 

Tables  4.1  to  4.6  and  Figures  4.1  to  4.3  are  the  results  of  the 
three  iterations  to  determine  the  optimal  policy.  The  control  policy, 
shown  at  the  top  of  each  figure,  is  the  one  used  during  the  forward 
iteration  resulting  in  the  states  shown.  The  adjoint  variables  on 
the  reverse  iteration  are  also  shown  at  the  bottom  of  each  figure.  The 
final  time  decreased  from  13.2  minutes  to  10.27  minutes  in  three 
iterations.  The  maximum  allowed  change  in  the  flow  rates  were  1  lb/mi n 
for  u-j  and  u~  and  10  lbs/min  for  the  recirculation  rate  u-,.  The 
maximum  allowed  change  in  the  steam  temperature  was  5°F.  The  minimum 
time  for  this  scenario  could  not  be  reduced  any  further  as  all  the 
variables  were  bounded. 

It  can  be  noticed  from  the  plots  that  the  final  control  policy 
is  bang-bang  and  that  the  switches  occur  at  the  point  constraints 
assumed  in  the  scenario.  This  is  indeed  a  fortunate  result  as  the 
control  policy  can  be  put  in  feedback  form  dependent  upon  whether  the 
states  are  below  or  at  their  steady  state  values.  This  is  not  a  typical 
result  but  is  due  to  the  assumed  control  scenario  which  is  specific  to 
this  type  of  problem.  The  final  control  policy  which  resulted  in  a 
minimum  time  can  thus  be  put  into  the  following  feedback  form. 
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Figure  4.1  Control,  State  and  Adjoint  Variables  for  Problem  1,  Iteration  1 


90 


u3  =  75.0 


lu. 


u1  =  u?   =  11 .0  II .0 


'u3  =  0     u2  =  0 


u4  -   245.0 


Ul  =  U2 


260r- 


-.40.0 


Figure  4.2  Control,  State  and  Adjoint  Variables  for  Problem  1,  Iteration  2 
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Figure  4.3  Control,  State  and  Adjoint  Variables  for  Problem  1,  Iteration  3 
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Intereffect  liquid  flow,  u2  =  u2  max  for  x-j  <  x-,  (first  effect  filling) 

Intereffect  liquid  flow,  u2  =  V21    for  Xi  >  x-|  (first  effect  filled) 

Feed  flow,  u-|  =  u-j  max  for  x3  <  x3  (second  effect  filling) 

Feed  flow,  u,  =  u2  +  VQ2  for  x3  >_  x3  (second  effect  filled' 

Recirculation  rate,      u3  =  u3  max  for  x-|  >  x-|  (recirculation  after 

first  effect  filled) 
Steam  temperature,      u4  =  u4  max  for  all  time  (input  steam  tempera- 
ture at  maximum) 
With  the  optimal  policy  as  shown  above  a  simulation  of  the 
state  equations  including  the  concentration  equations  was  done  and 
is  shown  in  Figure  4.4  while  Table  4.7  gives  the  simulated  values. 
The  product  flow  rate,  u5,  was  kept  shut  off  until  the  desired  concen- 
tration was  reached  in  the  first  effect.  While  concentration  was 
taking  place,  the  two  hold-ups  were  kept  steady  by  maintaining 

u2  =  V2]   and   U]  =  u2  +  VQ2 

After  the  desired  concentration  C,  was  achieved,  the  product  flow  rate 
was  fixed  at  u,-  =  u?C?/C,  to  keep  the  product  concentration  constant. 

In  the  simulation  shown  in  Figure  4.4,  it  was  assumed  that  the 
initial  and  final  concentrations  of  the  solute  were  very  small  so  as 
to  cause  a  negligible  boiling  point  elevation  and  also  so  that  the 
properties  of  water  were  not  altered  by  the  presence  of  the  solute. 
For  example,  a  3  percent  weight  solution  of  sodium  hydroxide  would  cause 
a  boiling  point  elevation  of  approximately  3°F  at  212°F  which  was 
considered  to  be  an  insignificant  elevation  for  this  problem. 
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Figure  4.4  Optimal  Simulation  including  Concentration  dynamics 
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IV. 3. 2  Problem  2.  Fixed  Feed  Rate 


The  minimum  start-up  time  was  again  determined  using  the  same 
algorithm  for  a  problem  similar  to  Problem  1  but  with  the  added 
restriction  that  the  feed  rate  was  fixed  at  a  nominal  value.  Such 
a  start-up  procedure  is  used  by  the  students  in  the  undergraduate 
unit  operations  laboratory.  It  was  desired  to  determine  the  effect 
of  this  policy.  (No  verification  was  made  that  industrial  practice 
uses  this  procedure,  but  it  would  be  surprising  if  many  start-up 
policies  did  not.)  The  nominal  value  chosen  was  2.9  Ibs/min 
(0.025  kg/s)  which  is  the  steady  state  flow  rate  used  for  the 
laboratory  experiment.  It  is  clear  that  the  minimum  time  must  be  at 
least  23  minutes  as  this  is  the  time  needed  just  to  fill  the  two 
effects.  The  product  flow  rate  from  the  first  effect  then  became 
an  active  control  variable  and  had  to  be  used  to  maintain  a  constant 
hold-up  in  the  first  effect.  The  control  scenario  was  modified  for 
this  problem  as  the  final  time  was  determined  by  the  hold-up  in  the 
second  effect  reaching  its  steady  state  value.  The  second  effect 
solution  started  boiling  at  an  earlier  time  and  so  was  not  used  as 
the  stopping  condition  as  in  Problem  1.  The  scenario  was  as  follows: 
Stage  A:  tQ   <  time  t  <  t, .  First  and  second  effects  are  being  filled. 
First  effect  is  being  heated. 
Control  variables: 

Recirculation  flow,  u3,  and  product  flow,  u5,  are  not  possible 
since  the  first  effect  hold-up  has  not  reached  its 
desired  value. 
Feed  to  first  effect,  u2,  and  steam  temperature,  u4,  are 
found  from  minimization  of  the  Hamiltonian. 
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Hamiltonian: 

Ha  =  A1U2  +  x  [u2(x4  "  X2J  +  Ql]  +  A3(WF  "  u2> 

The  second  effect  enthalpy  is  constant;  x,  =  0  or  x,  =  lv  resulting 

in  X«   =  0.  Time  t-, ,  signifying  the  end  of  Stage  A,  is  determined 

when  the  first  effect  is  filled,  i.e.,  when  x-,(t-,)  =  x, . 

Stage  B:  t,  <_  time  t  <  t~.  First  effect  is  filled  and  is  being 

heated.  Second  effect  is  being  filled. 

Control  variables: 

Feed  to  the  first  effect,  u?  =  0  to  maintain  a  constant 

hold-up.  Product  from  first  effect,  Ur  =  0.  Minimiza- 
tion of  the  Hamiltonian  determines  the  recirculation 
rate,  u~>  and  the  steam  temperature,  u«. 

Hamiltonian: 

Hb  =  -r-  Qi  +  *3wf 
xi 

The  first  effect  hold-up  and  the  second  effect  enthalpy  are 

unchanged;  x-,  =  x»  =  0;  and  A,  =  A,  =  0.  Time  tp»  signifying  the  end 

of  Stage  B,  is  determined  when  the  first  effect  liquid  starts  to  boil, 

i .e. ,  when  x2(t2)  =  x2. 

Stage  C:  t«  1  time  t  <  t,.  Second  effect  is  being  filled  and  being 

heated. 

Control  variables: 

The  feed  to  the  first  effect,  u2,   is  such  that  it  maintains 

a  constant  first  effect  hold-up,  u2  =  Vpi .  The  product 

rate  Ur  =  0.  The  recirculation  rate,  u3>  and  the 

steam  temperature,  u„,  are  determined  from  minimization 
of  the  Hamiltonian. 
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Hamiltonian: 

X 


H     =  ■£  CW12(h12  -  x2)  +  V21(x2  -  h\)  +  Q^  +  X3(WF  -  u2) 
xl 

XA 

+  ^  [WF(hF  -  x4)  +  Q2] 

Time  t,,  signifying  the  end  of  Stage  C,  is  determined  when  the  second 
effect  starts  to  boil,  i.e.,  when  x4(t3)  =  x4- 
Stage  D:  t3  <  time  t  <  tf .  Second  effect  is  being  filled. 
Control  variables: 

Same  as  in  Stage  C. 
Hamiltonian: 


H.  -  £  [W]2(h12  -  x2)  +  V21(x2  -  h])   +  Q^  +  X3(WF  -  u2  -  VQ2) 
xl 
First  effect  hold-up  and  second  effect  enthalpy  are  maintained  constant. 

Time  tf,  is  determined  when  the  second  effect  hold-up  reaches  its 
desired  value,  i.e.,  when  x3(tf)  =  x^. 

The  final  condition  on  the  Hamiltonian  yields  a  value  for 

A3(tf) 

X2(tf)  =  0,  x3(tf)x3(tf)  =  -1  determines  x3(tf) 

At  intermediate  points  continuity  of  the  Hamiltonian  is  used 
to  determine  the  unknown  multipliers.  At  time  t3, 

Hd(t3)  =  H  (t3)  determines  A4(t~) 
At  time  t, , 

Hb(t|)  =  Hfl(t[)  determines  X-j(t^) 

It  took  5  iterations  for  the  final  time  to  reach  its  minimum 
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value  of  24.53  minutes.  The  feed  rate  was  fixed  at  2.9  lbs/minute 
(0.025  kg/s).  The  control  variable  which  contributed  significantly 
to  this  reduction  in  the  final  time  was  the  intereffect  flow  rate  u2- 
The  state  variables  for  three  of  the  iterations  are  shown  in  Figure  4.5 
and  the  values  are  tabulated  in  Tables  4.8  to  4.10.  The  results  of 
the  iterations  were. 

Iteration  # 


u^ 

Final  Time 

1.7 

26.875 

2.0 

26.43 

2.2 

25.35 

2.4 

24.625 

2.9 

24.531 

1 

2 
3 

4 
5 

Again,  it  can  be  noticed  from  the  plots  that  the  switches  in 
the  control  variables  take  place  at  point  constraints,  i.e.,  the 
switches  can  be  directly  related  to  the  states  yielding  a  feedback 
control  policy.  The  control  policy  can  be  put  in  feedback  form  as 
follows: 

Intereffect  liquid  rate,  u2  =  Up  for  x]  <  x-j  (first  effect  filling) 
Intereffect  liquid  rate,  u2  =  V21  for  x]  >  x1  and  x3  <  x3  (first  effect 

filled,  second  effect  filling) 
Intereffect  liquid  rate,  u2  =  Wp  -  VQ2  for  x3  >  x3  (second  effect 

filled) 
Recirculation  rate,  u3  =  0  for  x-j  <  x^  (first  effect  filling) 
Recirculation  rate,  u3  =  u3  max  for  x-,  >  x]  (first  effect  filled) 
Steam  temperature,  u4  -  u4  max  for  x2  <  x2  (first  effect  heating) 
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ITERATION  5 

ITERATION  4 

ITERATION  1 


40. On 


235 


0.0    5.0     10.0     15.0     20.0     25.0     30.0 

TIME(MINUTES) 

Figure  4.5  State  Variables  for  Problem  2,  Iterations  1 ,  4  and  5 
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Steam  temperature,  u*  =  u*  .  for  x?  £  x2  (first  effect  boiling) 
Product  rate,  u^  =  0  for  x,  <  x.,  (second  effect  filling) 
Product  rate,  Ur  =  u~  -  V^-i  for  x~  >_  x,  (second  effect  filled) 

IV. 3. 3  Problem  3.  No  Bound  on  the  Second  Effect  Hold-up 

The  minimum  time  problem  was  run  again  for  the  conditions  as 
in  Problem  1  but  with  the  upper  bound  on  the  second  effect  hold-up 
removed.  The  second  effect  hold-up  H?  (state  variable  x~)  was  permitted 
to  vary  with  the  restriction  that  it  should  be  within  one  percent  of 
its  desired  steady  state  value  at  the  final  time.  Since  the  hold-up 
is  allowed  to  vary  the  feed  to  the  second  effect,  u, ,  need  not  be 
equal  to  the  flow  of  liquid  and  vapor  out  of  the  second  effect, 
Up  +  Vgp,  and  minimization  has  to  be  done  with  respect  to  u-,  and  u?  as 
opposed  to  Up  alone  as  in  problem  1.  The  control  scenario  is  different 
from  the  one  used  in  problem  1  as  we  now  assume  that  x4(tf)  =  x. 
and  x-Jt.)  =  x\  ±  0.01  x\  at  the  final  time.  In  other  words  two  of 
the  states  reach  steady  state  at  the  final  time  resulting  in  an 
iterative  solution  to  update  the  guess  on  one  of  the  adjoint  variables 
at  the  final  time.  The  scenario  is  as  follows: 

Stage  A:  tQ  <_  time  t  <  t-,  .  The  two  effects  are  being  filled  while 
the  first  effect  is  being  heated. 
Control  variables: 

The  first  effect  recirculation  flow,  u~,  and  product  flow,  Ur, 


are  not  possible  until  the  first  effect  hold-up  has 
reached  its  steady  state  value.  Thus  u3  =  u5  =  0. 
Minimization  of  the  Hamiltonian  is  with  respect  to  the 
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feed  to  the  second  effect,  u, ,  the  intereffect 
flow  rate,  u?,  and  the  steam  temperature,  u.. 
Hamiltonian: 

h 

Ha  =  Alu2  +  x~  ^u2^x4  "  x2^  +  °^   +  A3^ul  "  u2^ 

The  second  effect  enthalpy  is  unchanged,  x.  =  0,  x.  =  lv,  X.  =   0. 

Time  t, ,  signifying  the  end  of  Stage  A,  is  determined  when  the  first 

effect  hold-up  reaches  its  steady  state  value,  i.e.,  x-,(t-,)  =  x, . 

Stage  B:  t-,  <_  time  t  <  t?.  First  effect  is  being  heated  and  the 

second  effect  is  being  filled. 

Control  variables: 

The  feed  to  the  first  effect  is  stopped,  i.e.,  u2  =  0, 
to  maintain  the  first  effect  hold-up  at  x-, .  The 
product  rate  Ur  is  kept  off,  i.e.,  u,-  =  0.  Minimiza- 
tion of  the  Hamiltonian  is  with  respect  to  the 
second  effect  feed  rate,  u-, ,  the  first  effect  recircula- 


tion rate,  u.,,  and  the  steam  temperature,  u„, 


Hamiltonian: 


A2 
Hb  =  —  [u2(x4  -  x2)  +  Q-,]  +  X3(u-j  -  u2) 

xl 
The  first  effect  hold-up,  x-, ,  is  unchanged;  x,  =  0,  x,  =  x, , 
X-,  =  0.  The  second  effect  enthalpy  is  also  unchanged; 
x4  =  0,  x4  =  hp,  A4  =  0. 
Time  t?,  signifying  the  end  of  Stage  B  is  found  from  the  condition 
x2(t2)  =  x2,  that  is  when  the  enthalpy  of  the  solution  in  the  first 
effect  reaches  its  steady  state  value. 
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Stage  C:  t?  <_  time  t  <  tf.  The  second  effect  is  being  filled  and 

heated. 

Control  variables: 

The  feed  to  the  first  effect,  u?,  is  set  to  maintain  the 
first  effect  hold-up  constant,  u?  =  Vp-i  •  The  product 
rate,  Ur,  is  kept  at  zero.  Minimization  of  the 
Hamiltonian  is  with  respect  to  the  feed  to  the  second 
effect,  u, ,  the  first  effect  recirculation  rate,  u,, 
and  the  steam  temperature,  u». 

Hamiltonian: 

X2  v 

Hc  =  ~  [W12(h12  -  x2)  +  V2](x2  -  \\\)   +  Q^  + 

xl 

X3(u1  -  u2)  +  ■—  [u^hp  -  x4)  +  Q2] 

The  first  effect  hold-up  is  unchanged;  x-,  =  0,  x-,  =  x, , 
X-,  =  0. 
Time  tf,  the  final  time  is  determined  when  the  second  effect  solution 
starts  to  boil,  i.e.,  x,(tf)  =  x».  At  this  time  we  also  require 
that  the  second  effect  hold-up  reach  a  certain  value  x^(tf)  =  x^  i 
0.01*x3. 

This  change  in  the  scenario  requires  that  one  of  the  adjoint 
variables  at  the  final  time,  A3(tf),  be  guessed  initially  and  then 
updated  on  subsequent  iterations.  The  multiplier  A,(tf),  corresponding 
to  the  stopping  condition  x,(tf)  =  x»  (solution  starts  to  boil), 
is  found  from  the  final  value  of  the  Hamiltonian.  A3(tf)  is  updated 
on  successive  iterations  by  choosing  it  to  zero  the  gradient  of  L  with 
respect  to  it,  which  is  equivalent  to  driving  >L  -  x^(tf)  to  zero. 
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Tables  4.11  to  4.15  and  Figures  4.6  to  4.8  show  the  results 
of  three  of  the  iterations.  Of  interest  is  that  tf  is  decreased  from 
10.27  minutes  to  9.80  minutes.  The  initial  policy  was  taken  to  be 
the  final  policy  from  Problem  1.  On  the  reverse  iteration  A3(tf) 
(corresponding  to  the  second  effect  hold-up)  was  guessed  to  be  zero. 
This  resulted  in  a  control  policy  that  was  bang-bang  on  all  the 
variables  but  which  required  the  switching  time  (T )  on  the  second 
effect  feed  rate,  u-, ,  to  be  increased.  The  results  for  iterations  2 
through  7  are  summarized  below. 


Ul 

ul 

Iteration 

sw 
5.875 

t<Tsw 
1  ,max 

t>Tsw 
ul,min 

V 

x3-x3(tf) 
-2.9 

x3(tf) 

2 

10.375 

0.17 

3 

5.85 

Ut 

1  ,max 

ul  ,min 

10.25 

-2.0 

0.09 

4 

5.8 

ul ,max 

ul ,min 

10.13 

-0.34 

0.04 

5 

5.75 

ul ,max 

ul  ,min 

10.13 

-0.07 

0.01 

6 

5.7 

Ut 

1  ,max 

ul  ,min 

9.8 

0.0 

0.005 

7 

5.6 

ul  ,max 

ul  ,min 

9.75 

2.92 

The  second  effect  feed  flow  rate,  u, ,  was  switched  from  its 
maximum  value  to  its  minimum  value  at  time  T    Integration  of  the 

Sw 

state  equations  was  continued  until  the  second  effect  started  to  boil, 
i.e.,  when  x-(tf)  =  x\.  At  this  time  tf,  the  deviation  of  x3(tf) 
from  its  desired  value  x\  was  computed  and  a  correction  applied  to  A3(tf) 
as  shown  in  the  right-most  column.  Convergence  was  obtained  at  the 
switching  time  5.7  minutes  as  this  had  a  minimum  final  time  as  well  as 
a  minimum  deviation  of  the  second  effect  hold-up,  x3,  from  its  steady 
value,  x3.  Any  further  reduction  in  the  switching  time  did  result 
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Figure  4.6  Control,  State  and  Adjoint  Variables  for  Problem  3,  Iteration  1 
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Figure  4.7     Control,  State  and  Adjoint  Variables  for  Problem  3,   Iteration  3 
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in  a  shorter  final  time  but  at  the  expense  of  a  larger  deviation  of 

the  second  effect  hold-up,  x3,  from  the  steady  state  value,  x3. 

The  optimal  policy  for  this  problem  can  be  put  in  a  feedback 

form  by  getting  rid  of  the  dependence  of  the  feed  to  the  second  effect, 

u, ,  on  the  switching  time.  It  should  be  noted  that  the  switching  time 

for  u,  in  Problem  1  was  approximately  5.66  minutes  which  is  close  to 

that  for  Problem  2.  This  is  not  unusual  as  the  switch  from  u,  m,„ 

1  ,max 

to  u-,      ■     in  Problem  2  occurred  when  the  excess  hold-up  in  the  second 
1  ,rmn  r 

effect  was  just  sufficient  to  drain  out  by  the  time  the  second  effect 

liquid  started  to  boil.  As  the  vaporization  in  the  first  effect  is 

of  the  order  of  0.1  lbs  per  minute  (0.0008  kg/s),  the  second  effect 

should  be  overfilled  by  about  0.5  lbs  (0.23  kg),  above  the  final 

hold-up  of  32  lbs  (14.56  kg)  to  compensate  for  the  vaporization  in 

the  first  effect.  The  major  difference  between  the  optimal  policy 

of  Problem  3  and  that  of  Problem  1  is  that  the  second  effect  feed 

flow  rate  is  stopped  after  the  second  effect  fills  up  whereas  u-,  =  u2, 

the  intereffect  feed  flow  rate,  in  Problem  1  during  the  same  stage  of 

startup. 

Summary  of  the  control  policy: 

Feed  flow,  u0  =  u0  „   for  x-,  <  x-,  (first  effect  filling) 
c        Z, max     11  3 

Feed  flow,  u2  =  V21  for  x,  >_  x,  (first  effect  filled) 

Intereffect  liquid  flow,  u,  =  un  m,  for  x~  <  x~  (second  effect  filling' 

i    I  »max     o    -5 

Intereffect  liquid  flow,  u,  =  u-,   .  for  x3  >_  x3  (second  effect  filled) 
Recirculation  rate,  u,  =  0  for  x-,  <  x-,  (first  effect  filling) 
Recirculation  rate,  u,  =  uQ  m=v  for  x,  >  x\  (first  effect  filled) 
Steam  temperature,  u,  =  u,     for  all  time. 
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IV. 4  Experimental  Runs 

Start-up  runs  were  made  on  the  evaporator  as  a  final  test  of 
the  minimum  time  policy  with  a  bound  on  the  second  effect  hold-up 
(as  in  Problem  1,  Section  IV.3.1).  These  runs  consequently  served  as 
a  test  for  the  evaporator  model  as  well.  The  procedure  for  each  run 
was  as  follows: 

1 )  The  feed  to  the  second  effect  was  preheated  to  around 
180°F. 

2)  The  temperature  in  the  steam  chest  of  the  first  effect, 
u»,  was  kept  at  the  maximum  possible  value  throughout  the  run. 

3)  The  feed  rates  to  the  two  effects  were  at  their  maximum 
values  initially.  The  feed  to  the  first  effect  was  cut  off  as  soon 

as  the  first  effect  hold-up  reached  its  steady  state  value.  This  hold-up, 
X-, ,  was  then  put  on  analog  control  by  manipulation  of  the  feed  rate  to 
the  first  effect,  u?,  by  one  of  the  automatic  controllers. 

4)  The  recirculation  pump  was  started  and  the  recirculation 
flow  rate,  u~,  rapidly  built  up  to  its  maximum  value. 

5)  The  feed  rate  to  the  second  effect,  u-, ,  was  kept  at  its 
maximum  value  until  the  second  effect  hold-up,  x~,  reached  its  steady 
state  value.  This  hold-up  was  then  put  on  analog  control  by  manipulation 
of  the  feed  flow  rate  to  the  second  effect,  u-, . 

6)  The  flow  rates,  temperatures,  hold-ups  and  vacuum  pressure 
were  sampled  automatically  by  the  IBM  1070  unit  every   15  to  20  seconds. 
The  1070  unit  also  maintained  the  recirculation  rate  u3  constant  by 
moving  the  setpoint  of  the  controller  which  operated  on  the 
recirculation  valve  CV3  (see  Figure  2.2)  depending  upon  the  deviation 
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of  the  measured  flow  iu  from  the  desired  flow.  A  linear  Kalman  filter 
was  built  for  this  particular  flow  as  it  was  subject  to  a  lot  of 
pulsations  and  noise.  The  discrete  filter  was  formulated  as  follows: 

F  =  F  ,  +  v  ,  (4.26) 

n    n-1    n-1 

Y  =  F  +  w 
n    n    n 

where  F  is  the  flow  rate  u~  which  we  wish  to  maintain  constant.  Y 

is  the  measured  value  of  the  flow,  v  and  w  are  random  noise  sequences 

of  mean  zero  and  covariances  q  and  r,  respectively.  The  filter 

equations  were  obtained  from  Jazwinski  (1970). 

Given  F  ,  and  P  ,  ,  where  F  ,  is  the  estimate  of  F 
n/n     n/n       n/n 

and  P  ,  is  the  covariance  of  the  estimate,  we  can  predict  the  flow 
rate  and  covariance  to  the  next  sampling  time  by 


and 


Fn+l/n   Fn/n 


Pn+l/n  "  Pn/n  +  qn 


On  sampling  the  flow  rate  y  ,  at  the  next  sampling  time,  we  can  update 
the  covariance  and  estimate 

P2 
D      -  d         n+l/n 


n+l/n+1    n+l/n   P  ...  +  r  ., 
n+l/n    n+1 

F        -  F      +  Pn+1/n+1  (V  F      \ 

n+l/n+1    n+l/n    r.,   un+l  "  n+l/n; 


P 

"n+1 

This  was  the  recursive  scheme  used  to  obtain  an  estimate  of  the  inter- 
effect  flow  rate,  u3.  It  was  found  that  regardless  of  what  the  initial 
value  for  the  covariance  PQ  was  the  estimate  converged  rapidly—generally 
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within  three  or  four  sampling  points.  To  allow  for  this  convergence 
to  the  real  flow  rate,  control  was  not  done  at  the  first  few  sampling 
points. 

Figure  4.9  shows  the  behavior  of  the  filtered  flow  rate  u~ 
as  opposed  to  the  measured  flow  rate.  In  all  three  cases  shown  it  can 
be  seen  that  the  filtered  estimate  is  representative  of  the  actual 
flow  after  about  6  sampling  intervals  regardless  of  the  initial  estimate 
of  the  variance  of  the  estimate  p.  In  most  of  the  runs  an  initial 
estimate  of  unity  was  taken  for  p.  The  variances  of  the  process  and 
measurement  noise  were  also  taken  as  unity. 

The  results  of  complete  runs  through  start-up  are  shown  in 

Figures  4.10  to  4.16.  In  these  figures  the  experimental  minimum  time, 

t   ,  is  compared  with  either  the  theoretical  minimum  time  predicted 
exp      r 

by  the  model,  t  .,  using  idealized  controls  or  the  minimum  time 
predicted  by  the  model,  t  ,,  using  the  actual  controls  measured  in  the 
experiment.  The  initial  and  final  conditions  for  the  model  simulations 
were  the  same  as  those  of  the  experimental  runs. 

Figure  4.10  is  a  comparison  of  the  experimental  and  actual 
minimum  times  for  Run  CI.  The  data  and  simulated  values  are  listed 
in  Tables  4.16  and  4.18.  The  control  variable  values  used  for  the 
simulation  are  the  smoothed  values  obtained  in  the  experiment.  Note 
that  in  this  run  and  in  the  remaining  runs  the  hold-ups  measured  are 
not  quite  near  those  predicted.  Some  possible  explanations  for  this 
discrepancy  are  given  in  Chapter  V.  Figure  4.11  compares  the 
experimental  time  with  the  theoretical  minimum  time,  t  . ,  using 
idealized  controls--in  other  words,  the  control  variables  used  on  the 
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simulation  were  shut  off  and  on  precisely.  Tabulated  values  for  this 
simulation  can  be  found  in  Table  4.17. 

Figure  4.12  shows  the  results  of  experimental  run  C2  and  the 
simulation  using  the  actual  controls.  As  these  controls  were  more 
realistic  than  the  idealized  controls,  all  subsequent  simulations  used 
the  actual  controls.  Data  and  simulated  values  for  this  run  are 
listed  in  Tables  4.19  and  4.20. 

Run  C3  was  a  start-up  run  in  which  the  recirculation  rate 
was  held  at  around  100  lbs/mi n  as  opposed  to  80  lbs/mi n  for  runs  CI 
and  C2.  The  experimental  and  simulated  values  are  plotted  in  Figure  4.13 
and  listed  in  Tables  4.21  and  4.22.  A  similar  recirculation  rate  was 
used  in  run  C4  and  these  results  can  be  found  in  Ficiure  4.14  and  in 
Tables  4.23  and  4.24. 

Runs  C5  and  C6  were  different  from  runs  CI  to  C4  in  that  the 
recirculation  rate  was  maintained  at  65  Ibs/min.  Figures  4.15  and 
4.16  along  with  Tables  4.25  to  4.28  contain  the  experimental  and 
simulated  values  for  runs  C5  and  C6. 

It  can  be  seen  from  Figures  4.10  to  4.16  that  the  model  is 
fairly  representative  of  the  process.  There  is  about  a  15  percent 
discrepancy  in  the  prediction  of  the  final  time  by  the  model  to  the 
final  time  obtained  experimentally.  This  could  be  either  due  to 
deficiencies  in  the  model  or  in  the  experimental  setup.  The  next 
chapter  lists  some  of  these  deficiencies  and  a  few  suggestions  are 
made  to  overcome  them. 
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LIST  OF  IMPORTANT  VARIABLES  USED  IN  CHAPTER  IV 

State  variables: 

X-,  -  First  effect  hold-up  (H-,) 

x?  -  First  effect  solution  enthalpy  (h,) 

x~  -  Second  effect  hold-up  (H2) 

x»  -  Second  effect  solution  enthalpy  (h2) 

Xr  -  First  effect  solute  concentration  (C-.) 

Xr  -  Second  effect  solute  concentration  (C^) 

Control  variables: 

u-,  -  Feed  flow  rate  to  second  effect  (top) 

u?  -  Feed  flow  (intereffect)  rate  to  first  effect  (Wj2) 

u,  -  Recirculation  rate  in  first  effect  (W-j-j) 

u,  -  Steam  temperature  to  first  effect  (T$) 

ur  -  Product  rate  from  first  effect  (Wq-j) 

Other  variables: 

W,?  -  Flow  rate  of  solution  entering  first  effect 

V21  -  Vapor  rate  out  of  first  effect 

Vn?  -  Vapor  rate  out  of  second  effect 

h1?  -  Enthalpy  of  solution  entering  first  effect 

hV  -  Enthalpy  of  vapor  in  first  effect 

hi  -  Enthalpy  of  vapor  in  second  effect 

C,9   -  Solute  concentration  in  solution  entering  first 
u  effect 

Q,    -  Heat  transfer  rate  in  first  effect 
Q?    -  Heat  transfer  rate  in  second  effect 


CHAPTER  V 
COMMENTS  AND  RECOMMENDATIONS 

V.l  Model 


The  lumped  parameter  nature  of  the  model  causes  it  to  be 
somewhat  inaccurate  in  predicting  the  temperatures.  There  is  a 
temperature  gradient  in  both  effects  and  particularly  in  the  first 
effect  as  is  evidenced  by  the  temperature  gauges  at  the  top  and 
bottom  of  the  tubes.  Using  a  temperature  gradient  based  on  these 
two  temperatures  is  probably  better  than  using  a  bulk  temperature. 

The  heat  transfer  mechanisms  assumed,  particularly  the  two- 
phase  coefficients  and  the  natural  convection  coefficient  in  the 
second  effect,  are  doubtful.  Developed  correlations  in  the  litera- 
ture are  of  doubtful  accuracy.  The  method  for  calculating  an 
average  heat  transfer  coefficient  in  two-phase  flow  is  particularly 
suspect.  All  of  these  suffer  some  loss  of  credibility  during  highly 
transient  conditions.  It  was  hoped  that  the  parameter  estimates 
G-i  and  9?  would  account  for  these  inaccuracies,  but  experimental 
runs  C  described  in  Chapter  IV  indicate  that  the  first  effect  tem- 
perature in  particular  is  not  predicted  well  once  boiling  has  begun. 

A  better  model  fit  would  have  been  possible  by  taking  9,  and 
6p  to  be  time  varying  parameters,  and  these  could  have  been  optimally 
determined  in  the  same  fashion  as  the  optimal  control  policy  was 
determined  in  Chapter  IV.  The  time  varying  nature  of  the  parameters 


149 


150 

could  possibly  have  accounted  for  varying  heat  loss  rates  throughout 
the  experiment  as  well  as  for  varying  temperature  gradients  in  the 
effects.  All  this,  of  course,  would  be  at  the  expense  of  added  com- 
puter time. 

V.2  Experimental  Setup 

The  measurements  of  the  levels  in  the  two  effects  were  sub- 
ject to  a  few  inaccuracies.  Measurement  was  based  on  the  height 
of  water  column  in  the  sight  glasses.  The  actual  height  varied  de- 
pending upon  the  intensity  of  boiling  and  the  pulsating  flow,  par- 
ticularly in  the  first  effect.  This  drawback  could  have  been 
alleviated,  but  certainly  not  eliminated,  had  the  DP  cells  which 
were  used  on  the  hold-ups  had  an  adjustable  damping  device  to  "filter" 
the  readings. 

These  momentary  variations  in  the  indicated  height  caused 
improper  analog  control  of  the  hold-ups.  The  proportional  gain  on 
the  automatic  controllers  had  to  be  kept  high  enough  to  enable  the 
control  valves  to  respond  reasonably  fast  to  changes  in  the  hold-ups, 
but,  with  a  high  gain,  the  momentary  variations  caused  by  the  boiling 
in  the  two  effects  would  cause  the  controller  to  respond  to  the  noise. 
A  compromise  value  of  the  proportional  gain  was  used,  but  this  did 
not  completely  overcome  the  problem  of  sensitivity  to  the  momentary 
variations.  A  filter  on  the  hold-ups  and  direct  control  of  the  set- 
points  by  the  IBM  1070  is  a  possible  solution  to  the  problem.  This 
can  be  done  if  the  vapor  rates  from  the  two  effects  could  be  directly 
determined  or  accurately  estimated  without  a  great  deal  of  calculation 
as  this  increases  the  on-line  computer  costs. 
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The  recirculation  flow  rate  in  the  first  effect  was  one  which 
did  cause  a  number  of  problems.  Starting  the  recirculation  pump  was 
critical  as  a  certain  amount  of  liquid  head  on  the  suction  side  was 
necessary  to  avoid  causing  the  pump  to  suck  in  air.  Thus  the  hold-up 
had  to  be  maintained  above  a  minimum  value.  In  addition,  the  hold-up 
had  to  be  kept  below  a  certain  value  not  much  above  the  minimum  re- 
quired to  avoid  entrainment  of  the  liquid  with  the  vapor  going  to  the 
second  effect.  This  hold-up  was  critical  and  had  to  be  maintained 
throughout  the  run. 

This  problem  of  a  critical  hold-up  could  be  avoided  if  the 
vapor-liquid  separator  were  baffled  or  if  an  additional  well-designed 
separator  were  installed  in  addition  to  the  existing  one.  This  would 
also  ensure  no  entrainment  at  higher  recirculation  rates. 

Temperature  measurement  at  the  exit  line  of  the  second  effect 
is  not  a  proper  means  for  obtaining  a  measure  of  the  temperature  in 
the  second  effect.  This  measurement  is  accurate  only  if  there  is  a 
flow  out  of  the  second  effect.  The  thermocouple  should  project  directly 
into  the  bulk  of  the  second  effect  liquid  in  order  to  get  a  good  estimate 
of  the  temperature.  The  use  of  additional  thermocouples  at  the  top  of 
the  two  effects  may  help  to  give  a  better  temperature  average  in  the 
two  effects  and  could  account  for  temperature  gradients. 

The  minimum  time  runs  C  shown  in  Chapter  IV  required  the  bleed- 
ing of  the  DP  cells  while  the  run  was  in  progress.  This  led  to  in- 
accurate flow  measurements,  but  only  at  the  initial  time. 
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V.3  Theory 

The  theoretical  development  of  the  minimum  time  solution  pre- 
sented in  Chapter  IV  is  not  a  generalized  approach  as  it  draws  heavily 
on  the  concept  of  a  scenario  leading  to  particular  orderings  of  the 
point  constraints.  Proposing  an  optimal  scenario  is  intuitively 
obvious  in  some  instances.  Most  start-up  problems  in  chemical  engi- 
neering processes  can  be  solved  by  examining  relatively  few  scenarios 
to  determine  the  optimal  one.  This  approach  of  ruling  out  most  of 
the  less  likely  scenarios  makes  the  mathematical  problem  much  less 
formidable. 

In  the  case  of  the  double  effect  evaporator,  the  minimum  time 
solution  is  bang-bang  on  the  control  variables.  Although  a  bang-bang 
minimum  time  solution  is  the  rule  in  linear  systems,  it  is  not  neces- 
sarily true  in  the  case  of  nonlinear  systems  such  as  the  evaporator. 
The  switching  times  found  here  were  directly  related  to  the  states  thus 
obtaining  the  control  policy  in  feedback  form.  This  form  made  the 
optimal  policy  easy  to  implement  in  practice. 

V.4  Conclusions 

In  spite  of  all  the  drawbacks  mentioned  in  the  preceding 
sections,  it  can  be  concluded  that  the  dynamic  and  algebraic  equations 
of  Chapter  III  together  with  the  parameter  estimates  for  6-j  and  Q^ 
provided  a  reasonably  accurate  working  model  of  the  double  effect 
evaporator.  The  experimental  runs  of  Section  IV. 4  show  that  the  model 
is  accurate  to  within  15  percent  in  the  prediction  of  the  final  time. 
Moreover,  with  the  control  policies  for  minimum  time  start-up  in 
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feedback  form  as  obtained  in  Chapter  IV,  the  accuracy  of  the  model 
does  not  really  matter.  The  control  policy  was  solely  a  function 
of  the  states  and  not  explicitly  a  function  of  time,  so  that  the 
precise  state  trajectories  predicted  by  the  model  were  not  important 
from  the  viewpoint  of  implementation  of  the  control  policy. 

The  scenario  approach  adopted  in  Chapter  IV  for  obtaining 
the  optimal  control  policy  and  the  resulting  use  of  point  constraints 
not  only  simplified  the  mathematics  but  was  a  factor  in  obtaining  the 
control  policy  in  feedback  form.  This  approach  could  possibly  be 
used  with  success  in  other  start-up  problems  in  the  chemical  industry. 


APPENDIX  A 
HEAT  TRANSFER  EQUATIONS  AND  OTHER  RELATIONSHIPS 

In  addition  to  the  6  differential  equations  (3.1)  to  (3.6) 

derived  in  Chapter  III  and  the  4  algebraic  equations  (3.7)  to  (3.9) 

and  (3.11)  we  have  the  following  algebraic  equations  to  complete 
the  set  of  equations. 

A.I  Relation  Between  Temperatures  and  Enthalpies 
T,  =  f(h-,)  -  from  solution  data  (A.l) 

To  =  f(ho)  -  from  solution  data  (A. 2) 

hi  =  f(T-,)  -  from  steam  tables  (A. 3) 

ho  =  f(To)  -  from  steam  tables  (A. 4) 

T  =  f(h  )  -  from  steam  tables  (A. 5) 

j       =   f(h-,p)  -  from  solution  data  (A. 6) 

TF  =  f(hF)  -  from  solution  data  (A. 7) 

A. 2  Heat  Transfer  Equations—First  Effect 
The  mechanism  for  heat  transfer  is  as  presented  in  (Fair,  1960, 
1963a,  1963b),  (Hughmark,  1969)  and  (Tong,  1965).  It  is  assumed  that 
the  first  effect  can  be  divided  into  two  zones — the  sensible  heating  zone 
(subscripted  s)  and  the  vaporizing  zone  (subscripted  B).  Single-phase 
convective  heat  transfer  is  the  mechanism  in  the  sensible  heating  zone 
whereas  the  mechanism  in  the  vaporizing  zone  is  a  combination  of  two- 
phase  convective  and  nucleate  boiling  heat  transfer.  Slug  flow  is 
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the  predominant  flow  pattern  in  which  the  vapor  bubbles  coalesce 
into  slugs  of  vapor  rapidly  once  boiling  begins.  The  slugs  of  gas 
accelerate  more  rapidly  than  the  remaining  liquid  and  the  vapor  has 
a  higher  velocity  than  the  liquid  at  exit.  The  fractional  lengths  of 
the  sensible  heating  zone  and  the  vaporizing  zone  are  obtained  through 
a  consideration  of  the  pressure  drops  relative  to  the  total  pressure 
drop  in  each  tube.  The  total  heat  transfer  rate  is 

QT  =  Qls  +  Q1B  (A.8) 

A. 2.1  Sensible  Heating  Zone 

Qls  =  UlsAir^ATls  (A.9) 

where  the  mean  temperature  difference  AT,  is 

ATls  =  Ts  -  <T1  +  T12>'2  <A-10) 

A,  is  the  area  of  the  sensible  heating  zone. 


(Lls 


h 

«1.  "  fi?sAl  rf-  <T,  "  W  <A-11) 

and  . 

where  T  ,  and  T  0  are  the  inside  and  outside  wall  temperatures, 
wis     w2s 

respectively.  Also 

"is  ■  fiiswAi  rf-  <T„is  -  W  (AJ3) 

where  h,  ,  is  the  film  coefficient  due  to  the  tube  walls  and  fouling. 
I  sw 

The  outside  film  coefficient  is  obtained  from  the  Mussel t  equation 
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h°s  =   0.943 


hsf  plsf  g  X1s 
Llsulsf(Ts  ■  W 


_ 


hs  -   f(V 


0.25 


(A. 14: 


(A. 15) 


The  inside  film  coefficient  is  obtained  from  the  Dittus-Boelter  equation 
calculated  at  boiling  zone  conditions 


hn  =  ni  d77 


DiiGi  ] 
hb 


f c ibhb 


0.4 


(A. 16) 


For  Re  = 


ylb 


>  5000;   n,   ■  0.023,   n2  =  0. 


2000  <  Re  <  5000;  n1  =  0.0775,  n2  =  0.667 
Re  <  2000;  n-,  =  0.183,  n2  =  0.545 
The  overall  heat  transfer  coefficient  in  the  sensible  heating  zone  is 
given  by 


Jls 


(A. 17: 


-o   7TT  D 


10  J. 


"is   elhls  U       hlsw 


This  introduces  an  unknown  parameter  6-,  to  account  for  the  uncertainty 
of  the  constant  used  in  equation  (A. 16).  6]  =  0]a  when  no  boiling  is 
taking  place  and  0-,  =  elb  when  there  is  boiling.  Also,  when  there  is 
no  boiling  L-j  =  L,. 

The  length  of  the  sensible  heating  zone  is  given  by 

AT 


H  AP  Js 


'Is 


(£M£M£ 


(A. 18) 
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where   Tp    is  the  slope  of  the  vapor  pressure  curve  estimated  from 


s 
steam  tables 


§   |  ■  f(V  (A. 19) 


The  term   I  tt  I    is  obtained  from  a  heat  balance 


All        *DliNlth1s(TwlB-V 

AL 


W12C1B 


AP 
The  term  I  -rr-  I  is  estimated  from 


(A. 20) 


^]=-P1Bg  (A.21) 


The  relationship  between  the  properties  of  liquid  occurring  in 
equations  (A. 14),  (A. 16),  (A. 20)  and  (A.21)  at  the  corresponding 
temperatures ,  are 

Tlsf  =  Ts  "  °-75(Ts  -  W  *  klsf  plsf  ylsf 
Tlb-  (T^  +  T^/2-.k^,  ylb,  Clb 

Tl  *  P1B'  C1B 
A. 2. 2  Vaporizing  Zone 

The  length  is  found  from  the  remainder  of  the  total  length 

Hb  "  H  "  hs  <A-22> 

The  heat  transfer  rate  is 

QlB  =  UlBAl  if  AT1B  (A-23) 
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The  mean  temperature  difference  is 


Also 


and 


AT1B  -  Ts   "  Tl 


QlB  =   R?BA1   Ef  <Ts"W 


L 


QlB  -   filVl   Ef  <TwlB"V 


Q1B       hlBwAl    U,      ^Tw2B  "  TwlB^ 


(A. 24) 

(A. 25) 
(A. 26) 
(A. 27) 


The  outside  film  coefficient  is  again  obtained  from  the  Nusselt 
equation 


h!JB  =  0.943 


"  3    2 
klBf  plBf  g  Xls 
p 

-LlBylBf(Ts  "  Tw2b' 


0.25 


(A. 28) 


The  fraction  of  vapor  to  liquid  is 


12       W 


(A. 29; 


The  variable  x-jp  "is  introduced  for  simplifying  the  solution  procedure 


^12 


0.4  x 


12 


(A. 30) 


The  Lockhart-Martinell  i  parameter  X..,  which  represents  a  ratio  of 
kinetic  energies  of  liquid  and  vapor,  is 


'     1             w            ^ 

0.9 

r     v     1 

0.5 

r  y 

0.1 

X^   = 

1    -  x12 

PjB 

IB 

tt 

x12 

,  PTB  j 

v 

1  hB  i 

r    ,            .      .0.9,      V      >0.5,             ^0.1 

X'      = 

1  -  x12 

P1B 

^B 

tt 

x12 

P1B  j 

V 

,  yiR 

(A. 31) 


(A. 32) 
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For  the  steam-air  system  the  two-phase  convective  heat  transfer 
coefficient  is 

hnP  =  3-5hisl4J  (A>33) 

The  mass  velocity  is 

Gl  =  N-T-  (A-34) 

1  ,Nirit 

The  coefficients  a  used  in  the  calculation  of  the  average  inside 
film  coefficient  are  given  by 


a12=f(G,,Xtt)  (A.35) 

aj2=f(GrXit)  (A.36) 

A  mean  a  is 

a12  =  (a12  +  a\2)/2  (A. 37) 

The  film  coefficient  due  to  nucleate  boiling  is  given  by 

filnb=52<Tw1B-Tl>  (A-38) 

The  average  inside  film  coefficient  is  a  combination  of  the  nucleate 
boiling  coefficient  and  the  two-phase  convective  coefficient 


filB  ■  "l^lnb  +  filtp  (A-39; 


The  overall  heat  transfer  coefficient  is 


U1B  =  1— k (A.40) 

J_+       1  10  +  J_ 

filB       elbfilB  Dl1        ^IBw 


160 


where is  a  combined  resistance  for  the  tube  walls  and  fouling. 

hlBw 
The  temperature-property  relationships  are 


Tl  "*  M1B'  P1B'  ylB  for  1ic'uid 


T]Bf  =  Ts  -  0.75(Ts  -  Tw2B)  -  k1Bfi  P1Bf,  y1Bf  for  vapor 


A. 3  Heat  Transfer  Equations—Second  Effect 


A-,  =  f(T1) 


Q2  =  UgAgATg 
where  the  mean  temperature  difference  is  given  by 


also 


AT2  =  T1    -   (T2  +  TF)/2 

^2  =  ^VT1  "  Tw2> 
Q2  =  h2A2(Twl  -  T2) 

The  outside  film  coefficient  is  found  from  the  Nusselt  equation 


(A. 41) 
(A. 42) 

(A. 43) 

(A. 44) 
(A. 45) 
(A. 46) 


h2  =   0.943 


2ffD2f  ^     1 
L2y2f(Tl    "  V 


0.25 


'A. 47' 


The  inside  film  coefficient  is  found  from  the  natural    convection 


equation 


Gr2  = 


LJJP29¥TW1    "  V 


(A. 48) 
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Pr 


C2y2 


2    k, 


Y2  =  Gr2Pr2 


0.55  k0  ] 


f  0.13  kn  ) 


0  25      4  7 

Y2    for  10^  i  Y2  <  3.5  x  107 


Y2*33  for  Y  >  3.5  x  107 


The  overall  heat  transfer  coefficient  is 


(A. 49) 
(A. 50) 

(A. 51) 


U2  = 


1 


h°  e2h2  D2i   h2w 


(A. 52) 


The  correction  factor  0?  is  again  introduced  to  account 
for  uncertainties  in  the  inside  film  coefficient  equation  (A. 51) 
6?  =  9?  when  the  liquid  is  not  boiling  and  02  =  e2h  w^en  *'ie  l1a>uic' 
is  boiling.  The  properties  of  vapor  required  in  (A. 47)  are  found  at 
the  film  temperature 


T2f  =  Tl  "  °'75(T1  "  TW2)  ~*   k2f  p2f  y 


2f 


The  liquid  properties  in  (A. 48),  (A. 49)  and  (A. 51)  are  found  at  the 
mean  temperature 


T2  =  ^T2  +  Tw2^/2  *  k2'  p2'  M25  C2'  32 


APPENDIX  B 
LISTING  OF  COMPUTER  PROGRAMS 

This  appendix  contains  a  listing  of  the  computer  programs 
written  for  determining  the  minimum  start-up  time  control  policy  for 
the  double  effect  evaporator  described  in  Chapter  II.  These  programs 
utilize  the  algorithm  of  Section  IV. 2. 4.  Basically,  each  iteration 
consisted  of  two  integrations  and  for  simplicity,  two  programs  were 
run  separately  with  the  output  from  one  being  the  input  to  the  other. 

The  first  program  handled  the  forward  integration  of  the  state 
equations  using  either  an  initial  guess  on  the  control  policy  or  an 
updated  version  of  the  previous  control  policy  as  input.  Integration 
was  done  by  subroutine  HPCG  (IBM  Scientific  Subroutine  Package)  which 
is  based  on  a  predictor  corrector  method.  Subroutine  XDOT  supplied 
the  derivatives  of  the  state  variables  which  are  the  right  hand  sides 
of  the  differential  equations  (3.1)  to  (3.6).  The  results  of  the  for- 
ward integration  comprising  the  state  and  control  variables  against 
time  were  output  on  cards  and  on  the  line  printer.  The  input  to  the 
next  main  program  consisted  of  these  cards  and  the  times  at  which 
point  constraints  were  encountered  on  the  forward  pass  together  with 
derivatives  of  the  states  before  and  after  these  points  in  time. 

The  second  program  integrated  the  adjoint  equations  in  the 
reverse  direction  in  time.  The  Hamiltonian  was  also  minimized  at 
selected  points  to  obtain  the  optimal  control  policy.  The  final 
conditions  and  jump  conditions  on  the  adjoint  variables  were  evaluated 
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using  the  derivatives  of  the  states  from  the  forward  integration. 
Subroutine  HPC6  was  again  used  for  the  integration  of  the  adjoint 
equations  and  subroutine  LAMDOT  supplied  the  right  hand  sides  of 
the  adjoint  equations.  Subroutine  VA04A  (Powell,  1964)  was  the 
search  program  used  for  minimizing  the  hamiltonian. 

Evaluation  of  the  time  dependent  variables,  other  than  the 
state  and  control  variables,  was  done  by  subroutine  FUNCS.  The  film 
coefficients  were  estimated  in  the  following  subroutines  with  the 
associated  references. 

1)  The  overall  coefficient  was  evaluated  in  subroutine  HEAT  based 
on  the  method  presented  in  Fair  (1960)  and  Tong  (1965). 

2)  The  forced  convection  coefficient  was  evaluated  in  subroutine 
FCIF  using  the  Dittus-Boelter  equation  (Hughmark,  1969). 

3)  The  natural  convection  coefficient  was  evaluated  in  subroutine 
FCIN  using  the  natural  convection  equation  (McCabe  and  Smith,  1967). 

4)  The  natural  convection  coefficient  for  two  phase  flow  was  evaluated 
in  subroutine  FCINB  using  the  Rohsenow  equation  (Tong,  1965). 

5)  The  steam  side  coefficient  was  evaluated  in  subroutine  FCO  using 
the  Nusselt  equation  (McCabe  and  Smith,  1967). 
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C       MAIN  PROGRAM  FOR  INTEGRATION  OF  THE  STATE  EQUATIONS 

£*•««*•••#«*«*•««»«••••«•*««*•*****•**••*•*»•*»•••*•****•*•*••**• 

C      PURPCSE 

C         THIS  IS  THE  MAIN  PROGRAM  WHICH  FIRST  CALLS  SUBROUTINE 

C         INPUT  TO  INPUT  DATA.   IT  THEN  INITIALIZES  THE  STATE 

C         VECTCR  AND  CALLS  SUBROUTINE  HPCG  FROM  THE  IBM  SCIENTIFIC 

C         SUBROUTINE  PACKAGE  TO  INTEGRATE  THE  STATE  EQUATIONS. 

C         FINALLY,  IT  OUTPUTS  THE  RESULTS. 

C      SUBPROGRAMS  REQUIRED 

C         SUBROUTINE  INPUT 

C         SUBROUTINE  HPCG 

C 

DIMENSION  X(4),PRMT(  5),.DERY(  4),AUX(  16,  A) 

EXTERNAL    XCOT,OUT 

COMMCN/TIMEl/TTT(  120  ) 

CONNCN/STATES/XK 120 ) , X 2 { 120 ) , X3 ( 120),XM 120) 

C0MM0N/TIMER/TIME1120)f.IT4MEi  IMAX 

CC^yCN/STEADY/XHAT(4),XINITU),TF,HF,N0RDER,UlMAX,UlMIN, 
lU2MAX,U2MN,U3MAX,U3MlN,U4MAXrU4MlN 

COMMCN/CCNTRL/UH 120 ),U2( 120),U3( 120),U4(120) 
C 

C  THE    INPUT    SUBROUTINE     IS    CALLED    HERE.        IT    INITIALIZES    THE 

C  CCNMON    BLOCKS     CATAl,     CATA2    AND    STEADY. 

C 

CALL    INPUT 
C 

C  INITIALIZE  ARRAY  COUNTER  ITIME  AND  VARIABLES  TO  BE  INPUT 

C         TO  SUBROUTINE  HPCG 
C 

ITIME=1 

N=4 

PRMT(1)=0. 

PRMT(2)=2C. 

PRMT( 3)=1.0 

PRMT(4)=?.2 

CO  13  J=1,N 
10  X(J)  =  XIMTU) 

DERY(l)=3./8. 

DERY(2)=l./8. 

DERY(3)=3./8. 

DERY(h)=1./8. 

CALL  HPCG(PRMT,X,CERY,N,  I FLF , XDOT ,OUT, AUX  ) 
C 

C         OUTPUT  RESULTS  ON  LINE  PRINTER  AND  ON  CARDS 
C 

wRITt(6,200C) 

I M=ITIME-1 

WRITE(6,2010)(K,TTT{K),X1(K),X2(K),X3(K),X4(K),K=1,IM) 

WRITE(6,2020) 

WRITl(6,2u10)(K,TTT(K),U1{K),U2(K),U3(K),U4(K),K=1,IM) 
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WRITE  (7,  2040)  (TTT(K),X1(K),X2(K),X3(K),X4{K),UI(K),.U2(K), 
1U3(K)  ,U4(K) ,K=1,  IM) 
STOP 
2000  F0RMAT(lHlt20X, •STATE  VARIABLES  ON  FIRST  I TERATION •,//, 5X  , 

1 'STAGE  NC ,5X, 'TIME' ,5X, 'XI' , 10X, 'X2',10X, 'X3' ,10X, 'X4' ) 
2010  FORMAT (7X, I3,7X,F6.3,2X,F9.5,4X,F9.5,4X,F9.5,4X,F9.5) 
2020  F0RMAT(1H1,20X, 'CONTROL  VARIABLES  ON  FIRST  I TERATION ',// , 

15Xt 'STAGE  NO'  ,5X, 'TIME' , 5X,«U1', 10X, 'U2',10X, • U3 • , 10X , • U4« ) 
2040  F0RMAT(F6.2,2X,8F9.4 ) 
END 
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C*« 
C 

c* » 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


XCGT 

PURPOSE 

THIS  SUBROUTINE  SUPPLIES  VALUES  OF  THE  DERIVATIVES  TO 

SUBROUTINE  HPCG. 
USAGE 

CALL  XCCT(T,X,DERY) 

NOTEO   WILL  BE  CALLED  ONLY  BY  SUBROUTINE  HPCG 
PARAMETERS 

T      =  TIME,  INPUT 

X      =.  STATE  VECTOR,  INPUT 

CERY   =  VECTOR  OF  DERIVATIVES  OF  X  AT  TIME  T,  OUTPUT 
SUBPROGRAMS  REQUIREC 

SUBROUTINE  MAPI 

SUBROUTINE  FUNCS 

SUBROUTINE  XDCT ( T , X , C ERY  ) 

DIMENSION  X(l  ),DERY(  1),Y(8),  IPOSTU)  ,F(8),-XX(8),JPOSTU)  , 
1D(8) 

COMMON/INDEX/  I  (A) 

COMMCN/TIMER/TIME(120),ITIME, IMAX 

COMMCN/TEMPU/UTEMP(A) 

CGMMCN/STEADY/XHATU),XINITU),TF,HF,N0RDER,U1MAX,U1MIN, 
1U2MAX,U2MIN,U3MAX,U3WIN,U4MAX,U4MIN 

COMMCN/CCNTRL/UK  120  ),U2(  120  ),.U3( 120),U4( 120) 

THE  ELEMENTS  OF  VECTOR  I  INDICATE  WHETHER  THE  STATES 

HAVE  ARRIVED  AT  THEIR  STEADY  STATE  VALUES  XHAT. 

EXAMPLE   1(1)  =  1   IF  Xd).GE.XHAT(l) 

=  0   OTHERWISE 

DO  2  N=l,4 

I  (N)=0 
8  IF(X(N)  .LT.XHAT(N)  )GC  TO  2 

I (N)=l 
2  CCNTINUt 

THE  VECTOR  UTEMP  CONTAINS  VALUES  OF  THE  CONTROL 
VARIABLES  TO  BE  USED  IN  THE  PRESENT  TIME  STEP.   THE 
OPTIMAL  LAW  IS  UTILIZED  HERE  IN  FEEDBACK  FORM.   IF 
A  SIMPLE  CONTROL  LAW  IS  NOT  AVAILABLE,  THEN  UTEMP 
SHCULC  BE  THE  INTERPOLATED  VALUES  OF  THE  CONTROL 
VARIABLES  INPUT  IN  FEEDFORWARD  FORM 

UTEMP(1)=U1MAX 

IF(I(3).EC.1)UTEMP(1)=U1MIN 

UTEMP(2)=U2MAX 

UTEMP(3)=U3MAX 

IF(  I  (  1)  .EG.C)UTEMP(  3)=U3M,IN 

UTEMP(4)=U4MAX 
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SUBROUTINE  MAPI  IS  USED  TO  MAP  THE  STATE  AND  CONTRGL 

VARIABLES  INTO  A  VECTOR  D  WITH  THE  STATE  VARIABLES 

IN  THE  FIRST  FOUR  POSITIONS  IN  D  AND  THE  CONTROL 
VARIABLES  IN  THE  NEXT  FOUR 

40  DO  10  J=l,4 
10  IPOST(J)=J 

CALL  MAPl(XX,8,X,  IP0ST,.4,Y  ) 

DC  20  J=5,8 

JJ=J-4 
20  JPOST(JJ)=J 

CALL  NAP1(Y,8,UTEMP, JPOST,  4,  C) 

THE  VALUES  OF  THE  TIME-DEPENDENT  VARIABLES  (WHICH  ARE 
FUNCTIONS  OF  THE  STATE  AND  CONTROL  VARIABLES)  ARE  FOUND 
BY  CALLING  SUBROUTINE  FUNCS 

CALL  FUNCS(D,8,F,8) 

RE-MAP  THE  CONTROL  VARIABLES  WHICH  MAY  HAVE  BEEN 
ALTERED  BY  SUBROUTINE  FUNCS 

DO  25  J=l,4 
25  UTEMPl J)=C( J  +  4) 
W12=F(1) 
V21=F(2) 
V02=F(3) 
Hl2=F(4) 
G1=F(5) 
G2=F(6) 
H1V=F17) 
H2V=F(8) 

THE  RIGHT  HAND  SICES  OF  THE  STATE  EQUATIONS  COMPRISE 
THE  DERIVATIVES 

DERYU)=W12-V21-UTEMP(3) 

IF(Xll) .LT..1E-70JG0  TO  50 

DERY(2)=(W12»(H12-X( 2 ) ) +V21* ( X( 2 )-HlV ) +Q 1 ) /  X< 1) 

GC  TC  60 
50  DERY(2)=0. 

GO  TG  70 
60  DERY(3)=UTEMP(1)-UTEMP(2)-V02 
70  IF(X(3)  .LT..DGO  TO  80 

DERY  (4)=  (UTEMPl  1)*<HF-X(  4)  ) +Q2+V02* I  X ( 4)-H2V )  )/  X  (3) 

IF(I (2).NE.1)CERY(4)=0. 

RETURN 
80  DERY(4)=0. 

RETURN 

END 


C  OUT 

£••««»«»*«*«*••**««***»«*•««••••**«*••«••*••*»*•*•••*••»»*•***••• 

C      PURPCSE 

C         THIS  SUBROUTINE  OUTPUTS  VALUES  OF  THE  STATE  VARIABLES 

C         ANC  DERIVATIVES  AT  TIME  T 

C      USAGE 

C         CALL  OUT(TT,X,DERY, I  HL F,ND IM,PRMT ) 

C         NCTEO   THIS  SUBROUTINE  IS  CALLED  ONLY  BY  HPCG 

C      PARAMETERS 

C         ALL  ARE  INPUTS.   REFER  TO  HPCG  FOR  DESCRIPTION 

C 

SUBROUTINE  OUT { TT , X, CERY ,  IHL F rND IM , PRMT ) 

DIMENSION  X(l ) ,DERY( 1),PRMT( 1) 

C0MMCK/STATES/X1(12C),X2(120),-X3(120),X4{120) 

COMMCN/CCNTRL/UK  120),U2(  120),.U3(  120),U4(120) 

CQMMCN/TIMER/TIME(120), I  TIME, IMAX 

COMMOK/TIMEl/TTTt 120) 

C0MMCN/STEADY/XHAT(4),XINIT(  4  )  , T F, HF , NORDER , U1MAX, U1M IN , 
1U2MAX,U2MIN,U3MAX,U3MIN,U4MAX,.U4MIN 

COMMON/TiZMPU/UTEMPJA ) 

COMMON/INDEX/ 1(4) 

I  1  =  I T IME 


THE  ERROR  PARAMETER  PRMT<4)  IS  RE-ESTIMATED  DEPENDING 
ON  THE  MAGNITUDE  CF  X  AND  THE  WEIGHTING 

PRMTU)  =  (0.375»ABS(X(  1)  ) +0. 125»ABS ( X { 2 ) ) +0. 3  75»AB S ( X ( 3 )  ) 
1+0.125*ABS(XU) ) )*0.01 

A  POINT  T(I)  IS  SAVED  IN  THE  OUTPUT  ARRAYS  (COMMON  BLOCK 
CCISTRL,  STATES  AND  TIMED  ONLY  IF  T(I)-T(I-2)  IS  GT.0.25 

MIKUTES 

IF( I  I  .LT.3)G0  TO  10 

J=II-1 

JJ=I  1-2 

IF(TT-TTTUJ)  -0.25)5,  10,  10 
5  K=J 

GC  TC  20 
10  K=  I  I 

11=11+1 

IF(  II.LT.12DG0  TO  20 

^RITE(6,2010)  I  I 

STOP 
20  X1(K)=X( 1) 

X2(K)=X(2) 

X3(K)=X(3) 

XMK)=.XU) 

Ul (K)=UTtMP(l  ) 

U2(K)=UTEMP(2  ) 


169 


2010 
2020 


U3(K)=UTEMP(3  ) 

U4<K)=UTEMP(4) 

TTT(K)=TT 

I TI ME=I I 

IFIK.EG. JJRETURN 

WRITE  (6,  2020  )TT,  (X(  J  ),DERY(  J  ),.J  =  1,NDIM)  , { PRMT < K ) , K=l , 5 ) 

I F ( I U).EC.1)PRMT(5)=1. 

RETURN 

F0RMATI5X, 'VECTOR  OV ERFLOWt 0 IMENS ION  =',I4) 

FORMATdOX, 'OUTPUT  FROM     OUT  AT  T I ME= • , F 10 . 5 , // , 5X , 
1*  STATES'  ,10X»  '  CERIVATIVES*f.4(/f.3X»FLC.5v  8X,E15.7)  ,/,lX,'PRM 
2T  PARAMETERS'  ,5(3X,F10. 5)  ) 

END 
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C»*» 

c 

c»»* 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


MAIN  PROGRAM  FOR  INTEGRATION  OF  THE  ADJOINT  EQUATIONS 

PURPOSE 

THIS  PROGRAM  CALLS  THE  INPUT  ROUTINE  TO  INITIALIZE 
COMMON  BLOCKS  CATA1.DATA2  AND  STEADY.   IT  THEN 
CALLS  CN  SUBROUTINE  HPCG  TO  INTEGRATE  THE  ADJOINT 

ATIONS.   WHEN  A  POINT  CONSTRAINT  IS  ENCOUNTERED  ON 
THE  REVERSE  INTEGRATION  CONTROL  RETURNS  TO  THIS 
PROGRAM  WHICH  DETERMINES  THE  VALUES  OF  THE  ADJOINT 
VARIABLES 

SUBPRCGRAMS  REQUIRED 
SUBROUTINE  SEARCH 
SUBROUTINE  CONVAR 

REAL  LAMDA(A) 

CI  MENS  I  CN  PRMT(5  )  ,  CL  AMDA  (  4  ) ,  AUX  (  16 , 4 ) „DER Y K ) 

EXTERNAL  L AMDCT, L AMCLT 

CCMMCN/PCINTS/TFtT2tTl,F3TFM,F4TFMtF2TlPtF2TlM,F3TlP,F3Tl*S 
1F1T1M 

COMMCN/STATES/XK  12C  ) ,  X  2  (  120  )  ,.X3  (  120  )  ,  X4  {  120  ) 

COMMCN/STEADY/XHATU),XINIT(  4  )  ,  TF  ,  HF  ,  NORDER  ,  U1MAX  ,  UlMIN  , 
lU2MAX,U2MIN,U3MAXtU3MIN,U4MAX,U4MIN 

COMMCN/CCNTRL/UH  120  )tU2<  120),.U3(120)  ,UM120) 

C0MMCN/TlMER/TlME(12o),ITIME,IMAX„IC0UNT,TIM,IG0 

COMMON/ C\EW/U1NEW< 120), U2NEW( 120 ) *U3NEW ( 120 ) ,U4NEW(120) 

CCCMCN/DITIME/NITIME 

READ  IN  VALUES  OF  THE  STATES  AND  CONTROLS  ON  THE  FORWARD 
INTEGRATION  INTO  COMMON  BLOCKS  STATES ,CONTRL  AND  TIMER. 
ALSO  THE  TIMES  AT  WHICH  POINT  CONSTRAINTS  ARE  ENCOUNTERE 

D 
ANC  THE  DERIVATIVES  BEFORE  AND  AFTER  THE  POINT  CONSTRAIN 

TS 

CALL  INPUT 

READ(5,1000) IMAX 

REAC(5,1010)(TlME(K),Xl(K)tX2(K),X3(K)fX4(K)fUl(K),U2(K), 
1U3(K)  ,U4 (K)  ,K  =  1, IMAX  ) 

READ(5,10  20)F3TFM,F4TFM,T2,T1,F2T1P,F3T1P,F1T1M,F2T1P, 
lF2TlMtF3riM,TF 

WRITE(6,2030) 

WRITE (6, 2040) (K,TIME(K),X1(KUX2(K),X3<K),X4(K),U1(K), 
1U2(K),U3(K),U4(K),K=1,IMAX) 

THE  FOLLOWING  PROGRAM  IS  FOR  A  SCENARIO  AS  IN  PROBLEM  3, 
CHAPTER  A.   GUESS  FINAL  CONDITION  ON  LAMDA13)  =  0.008 
FINAL  CONDITIONS  FOR  THE  ADJOINT  VARIABLES.   NOTE  THAT 
LAMCAU)  IS  FOUND  FROM  THE  FINAL  CONDITION  ON  THE 
HAMILTCNIAN 
IGCTl.  =  1 
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LAMDM1  )  =     • 

LAMDA(2)=0. 

LAMCA(3)=C008 

LAMDA(4)=(-1.+F3TFM*LAMDA( 3)  )/F4TFM 

ISTARTMMAX 

IEND=1 

PRMT(1)=TIME(  IMAX  ) 

PRMT(2)=0. 

NITIME=3 

PRMT(3)=-0.25 

PRMT(4)=0.0025 

N  =  4 

ITIME=ISTART 
10    CONTINUE 

DO    20    J=l,4 
20    OLANCM  J)=0.25 

IGO=IGOTC 

CALL    HPCG(PRMT,LAMDA,DLAMCA,Nt IHLF ,L AMDOT ,LAMOLT , AUX ) 

IF(IHLF.GT.10)STOP 
C 
C  CUTPUT    RESULTS 

WRITE(6,200C) 

WRITE(6,2010)(J*TIME(J),X1(J),.X2{J),X3U),X4(J),U1(J), 
1U2( J) ,U3{ J) ,U4( J ), J=  I  ST ART,  I  END) 

IGCTC=IGCTC+1 
C 

C         START  NEXT  INTEGRATION  AT  TIME  WHEN  PREVIOUS  STAGE  ENDED 
C 

PRMT( 1)=TIM 

GO  TC  (70,40*50,70),  IGOTO 
C 

C         CONDITION  CN  LAMDM4)  AT  TIME  T2 
C 

40  CONTINUE 

LAMDA(4)=G. 

GO  TC  10 
C 

C         CONDITION  ON  LAMCA(l)  AT  TIME  Tl 
C 


50  CONTINUE 

LAMDA(1)=( (F2T1P-F2T1M)*LAMDA(2)+(F3T1P-F3T1M)* 
1LAMDA(3)  1/F1T1M 
GO  TC  13 
70  CONTINUE 
STOP 
1000  F0RMATU3) 
1010  F0RMAT(F6.0,2X,8F9.0  ) 
1020  F0RMATI11F7.0) 

10  FORMAT! 1H1.20X, 'STATE  AND  CONTROL  VARIABLES', 

1         //,2X,'STAGE',2X,'TIME',5X,  'XISIOX,  •  X2  •  ,  10X,  •  X3  % 
210X, 'X4'  ,10X, «U1'  , 10X, 'U2', 10X, «U3', 10X,  'U4'  ) 
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2010  FORMAT (3X,I3t2X, F6.3,2X,F6.3,5X,F7.3,5X,F7.3,5X,F7.3,5Xf 

lF7.3,5X,F7.3,5X,F7.3,5XtF7.3) 
2030  FORMAT! 1H1,30X,« STATE  AND  CONTROL  VAR I A8LES  •  ,  //♦ 2 X, • STAG ■•  , 

15X,  »TIME«  t10X,»Xl«,lCX, •X2% 10X,  'X3« t10X, 'X4',10X, 'Ul1 ,10X, 

2'02'  .lOX.'US'  ,10X,«U4'  ) 
2C40  F0RMAT(4X,I2,5X,F6.2,5X,F7.2,5X,F7.2,5X,F7.2,5X,F7.2,5X, 

1F7.2,5X,F7.2,5X,F7.2»5X»F7.2) 
END 
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C  LAMOUT 

C      PURPOSE 

OUTPUTS  VALUES  OF  LAMCA  AND  THE  DERIVATIVES  AT  DESIRED 
TINES 
USAGE 

CALL  LAMOUT (TT , LAMCA, DERY , IHLF , ND IM , PRMT ) 
TC  BE  CALLEC  ONLY  BY  SUBROUTINE  HPCG 
PARAMETERS 

REFER  TO  SUBROUTINE  HPCG  FOR  PARAMETER  DESCRIPTION 

SUBRCUTINt  LAMOUT  (TT,  LAMCA,  DERY,.  IHLF,  NDIM,  PRMT) 
REAL  LAMDA(4),L1(4),L2(4),L3(4),L4(4),L(4) 
DIMENSION  CERY(l)  ,PRMT(  1),T(  4) 

C0MMCN/P0INTS/TF,T2,T1,F3TFM,F4TFM,F2T1P,F2T1M,F3T1P,F3T1M, 
1F1T1M 
CGMMCN/DITIME/NITIME 

COMMCN/TIMER/TIMEC 12C), I  TIME, IMAX, ICOUNT , TIM , I GO 
COMMCK/INCEX/ 1(4) 
COMMON/ CNEW/U1NEW( 120  )  , U2NEW ( 120  )  , U3NEW{ 120 ) ,U4NfcW(120) 

CHECK  IF  POINT  CONSTRAINT  IS  ENCOUNTERED.   IF  IT  IS 
ENCOUNTERED  TERMINATE  INT  GRATION  AND  RETURN  TO  MAIN 
PROGRAM  FOR  NEW  VALUES  OF  ADJOINT  VARIABLES 

IGCTC=IGO 

GO  TC  (60,70,80) , IGOTO 
60  CONTINUE 

IF(TT.LT.T2)PRMT(5)=1. 

GO    TO    80 
70    CONTIISUc 

IF(TT.LT.T1)PRMT( 5>  =  1. 
80  CONTINUE 

PRMT(4)=0.0025*{ ABS(LAMDA( 1) )+ABS(LAMDA( 2) ) +AB S t LAMDA ( 3 ) )+ 
1ABS(LAMDA(4) ) ) 

IF(PRMT(4).LT..0009)PRMT( 4)=0.C0C9 

TIM=TT 

T(  ICCUNT)=TT 

L1(ICCUNT)=LAMCA( 1) 

L2(  ICCUNT)  =  LAMDA( 2) 

L3(ICCUNT)=LAMCA( 3) 

L4(ICCUNT)=LAMDA(4) 

IF( ICCUNT.GT.3JG0  TO  10 

1CCUNT=ICCUNT+1 

RETURN 
10  TOUT=TIME( ITIME) 

IFITCUT. GE.T(4). AND. TOUT. LE. T( 1) )G0  TO  30 

DO  20  J=l,3 

T( J)=T( J+l) 

Ll( J)=L1( J+l) 
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L2( J)=L2( J+l) 
L3( J)=L3( J+l) 

20  L4( J)=L4( J+l) 
RETURN 

30  CALL  SEARCMT,  LI, TOUT, A,  3,L(  1)  ) 
CALL  SEARCH{T,L2,T0UT,4, 3,H 2)  ) 
CALL  SEARCMT,L3,T0UT,4,  3,L(  3)  ) 
CALL  SEA3Ch(TrL4,T0UT,4, 3,L( 4)  ) 

CALL  CCNVAR  ONLY  AT  TIMES  AT 
:         ON  THE  FORWARD  INTEGRATION 

CALL  CCNVAR(TCUT,L) 
WRITE(6,2Q40)T0UT,(L(J), J=l,4) 
WRIT £(7, 2040) TOUT, ( L  (  J  ) ,  J  =  I,  4  ) 

2040  F0RMAT(F7.3,3X,4F10.4) 
2C30  FCRf"AT(  I  2  ,  F8  .  3  ,  9  F8  .  3  ) 

ITIME=ITI^E-NITIME 

GO  TC  10 

END 


WHICH  VARIABLES  WERE  STORED 
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£»«««•«**»*»#«»****«#***«* 


£»«««««***«•«*««•#»«*•**»* 


LAMDOT 


PURPOSE 
PRCVI 
SUGRC 

USAGE 
CALL 

PARAMETE 
TT 

LAMDA 
DERY 

SUBPRCGR 
SUBRG 
SUERG 
SUBRG 
SUBRO 

SUBROUTI 
REAL  LAM 
CIMENSIC 


DES  DERIVATIVES  OF  THE  ADJOINT  VARIABLES  TO 
UTINE  HPCG 


LAMCOT(TT,LA 
NGTEO  CALLED 
*S 

-  TIME,  MI 

-  INPUT  VE 

-  OUTPUT  V 
AMS  REQUIRED 
UTINE  SEARCh 
UTINE  FUNCS 
UTINE  MAPI 
UTINE  PDERIV 

NE  LAMDOT(TT 

DA14) 

N  CERY(l) ,XX 


VCA,CERY  ) 
ONLY  BY  SUBROUTINE  HPCG 

NLTES,     INPUT 

CTOR    OF    ADJOINT    VARIABLES 

ECTOR    OF    DERIVATIVES 


lIP0SITI2)f IMEM(4) 
COMMCN/STATES/XK 120 
CCMM0N/CELT/DELX1,0E 
C0MMCI\/TIMER/TIME( 12 
COMMON/INCEX/  IU) 
CGNMCN/TEMP/X(8) 
CCCMCN/CCNTRL/UK 120 
COMMON/STEADY/XHATU 

1U2MAX,U2MIN,U3MAX,U3 
WRITE(6,2060)TT,LAMC 
2G60    FCRMATdOX, 'LAMDOT. 


114) 


INTERPOLATE  FOR  V 
VARIABLES  AND  SET 
BE  EXPRESSEC  IN  F 


CALL 

SEARCKTIMEtXl, 

TT 

CALL 

SEARCH(TIME,X2, 

TT 

CALL 

SEARCH(TIME,X3, 

TT 

CALL 

SEARCH<TlME»X<t, 

TT 

CALL 

SEARCMTIME,U2, 

TT 

DC    1 

J  =  5,8 

IF(X( J).LT.O.  )X( J  )  =  G 

. 

IF(X13).LE.Q.  )X(3)=G 

.07 

IF(X(  D.LE.C.  )X(  1)=0 

.001 

DC    2 

N=l,4 

,LAMCA,DERY) 
(2),FF(2,2),DFDX(2),W(2),F(8),D(8),Y(2) 


),X2«120),-X3(  120)  ,X4(120) 

LX2,DELX3,DELX4,0ELU1,DELU2,DELU3,DELU4 
C  ),.ITIME,  IMAX,  ICOUNT 


)tU2(  120)t-U3(  120)fU4(  120) 
),XINIT<  A),TFf.HF,NORDER,UlMAX,UlMlN, 
MINtUAMAX,U4MIN 
A, IMAX 

TT=«  ,  F10.3,  *LAMDAS=.'  ,4(2X«FL0.6)t2X,fIM 

AX=«  , 


ALUES  OF  THE  CONTROL  AND  STATE 

THOSE  CONTROL  VARIABLES  WHICH  CAN 
EEDBACK  FORM 


,  IMAX,N0RDER,X( 1)  ) 
,  IMAX,N0RDER,X(2) ) 
,  IMAX,N0RDER,X(3) ) 
,  IMAXf-NORDERt  X(4)  ) 
,  IMAX,N0RDER,X(6)  ) 
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IIN)=0 

IF(X(M .LT.XHAT(N) )GC  TO  2 

I IN)=1 

CCNTINUE 

X{8)=U4MAX 

X(5)=U1MAX 

I F  < I (3).EC.1)X(5)=U1MIN 

X(7)=U3MAX 

I F ( I (1) .Et.0)X(7)=U3NIN 

FEED=X(5) 

CALL  FUNCS(X,8,F,8) 

wl2=F(l) 

V21=F(2) 

VC2=F(3) 

H12=F(4) 

G1=F(5) 

G2=F(6) 

H1V=F(7) 

H2V=F(8) 

IGCTC=I( 1)+2*I(2)+4*I(3 )+8»I (4)+l 

GC  TG  (10, 170,50, 170,50, 170, 5C, 170, 50,50, 50, 150, 50, 


1IGCTC 

KC  STATE  VARIABLES  AT  STEADY  STATE 
DETERMINE  RIGHT  HAND  SICE  OF  CtLAM 


50,50,12 
5), 


MDAD/DT 


10  DERY(4)=0. 

13  XX(1)=X(1) 
KCUNT=0 

CALL  PDERlV(XX,FF,2,2,KOUNT,DELXl,DFDX,id) 
IF(KCUNT)5000, 16,  12 
Y( 1)^XX(KCUNT  ) 


11 
12 


I  POSIT (1 )  =  1 

CALL  l"APl(X,8»Yt  IPOSIT,  1,  D) 

CALL  FUKCSIC, 8,F,9) 

FF(KCUNT,1)=F(5) 


16 


FF(KCUNT,1)=F(5) 
FF(KCLNT,2)=F(6) 
GC  TC  11 
CC1DX1=CFCX( 1 ) 
DG2DX1=CFCX(2 ) 


1*(X(2)- 
10X1) /X( 

1)  + 
/X<3) 
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109    XX(1)=X(2) 
KCUNT=C 

106  CALL  PDERIV(XX,FF,2, 2 , K0UNTrD ELX2 ,DFDX , W ) 
IF(K0UNT)5000,108, 107 

107  Y(1)=XX(KCUNT) 
IPCS1T(D=2 

CALL  fAPKX,  8,Y,  IPOSIT,  1,  D) 
CALL  FUNCS(C,8,F,8) 
FF(K0UNT,1)=F(5) 
FF(K0UNT,2)=F(6> 
GC  TC  106 

108  DG1DX2=CFCX(1) 

3000  FORMAT(15X,'LAMDA(2)  =  'iF10.5,5X,.  •0Q1DX2=,,E15.7) 
DG2DX2=CFCX(2 ) 
IF(V21.LE.0.)CQ2DX2=0. 
DH1CX2=0. 

IF(W12.GT.C.)CH1DX2=X(7)/W12 
DENR=1078.-0.6*D(  2) 

DV2DX2=(CENR*CG2DX2+C.6*Q2)/(DENR*DENR) 
DHVDX2=0.4 

DV0DX2=CQ2CX2/(1078.-0.6*X(4 ) ) 
IFIVC2.LE.0.) CVODX2=0. 

DERY(2)=-LAMDA(2)»(Wl2*( DH1DX2-1. ) +DV2DX2* ( X ( 2 )-Hl V) 
1V21*(1.-CFVDX2)  +  0G1DX2)/X( 1 ) -LAMDA U ) * ( DQ2DX2  +  DVODX2 
2H2V))/X(3)+  (LAMDAl 1)»0V2DX2+L 


*(X(4)- 

AMOA(3) 


3DVCDX2) 

DETERMINE  RIGHT  HAND  SIDE  OF  D(LAMDA4)/DT 

155  CONTINUE 

XXil)=X(4) 

K0UNT=0 
160  CALL  PCERIV(XX,FF,2,  2..K0UNT,  DELX4,DFDX,W  ) 

IF(KCUNT)5000,180,165 
165  Y(1)=XX(KCUNT) 

IP0SIT(1)=4 

CALL  MAPl(X,8f.Yf  IPOSIT, ltD) 

CALL  FUNCS(C,8,F, 8) 

FF(KCUNT,1)=F(5) 

FF(KGUNT,2)=F(6) 

GO  TC  160 
180  DG1DX4=CFCX(1 ) 

DG2DX4=CFCX(2) 

IFU21.LE.0. )CQ2DX4=0. 

DH1CX4=0. 

IF(W12.GT.0. )CH1DX4=X(6)/W12 

DV2DX4  =  CG2DXW(  1078.-0.  6*X(  2)  ) 

DENR=1078.-0.6*D(  4) 

DV0DX4=( CENR*( CQ2CX4-FEEC ) +0 . 6* ( Q2-F EED* ( X( 4 )-HF ) ) ) / 


(DENR+D 
ENR) 
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IF(VC2.LE.0.)CVODX4=C. 

CHVDX4=0.4 

DERY{4)  =  -LAI"DA(2)»(W12*DHlCX4*(X(2)-HlV)»DV2DX4-«-.DQlDX4)/X(l 

)- 
lLA^OAK) »(-X(5)+DC2CX4*(X(4)-H2V )  »DV0DX4  +  V02» ( 1.-DHVDX4) )/X 

(3)  + 
2  (LAMDA(  l)*CV2DXA  +  LAMDA(3)*DVODX4) 

IF(  I  (M.£C.1)CERYK)=C. 

IF(  I  (2) .EQ.O)DERYi4)=0. 

GC  TC  112 

XI,  X2  AND  X4  AT  STEACY  STATE 

150  CCNTINUE 

DERY(1)=0. 
OERY(A)  . 
GC  TC  109 

XI  ANC  X2  AT  STEACY  STATE 

170  CONTINUE 

DERY( 1)=0. 
GC  TC  1G9 

DETERMINE  RIGHT  HAND  SICE  OF  C(LAMDA3)/DT 

112  XX(1)=X(3) 
KOUNT=0 

113  CALL  PCERIV(XX,FF,.2, 2 , KOUNT, CELX 3 , DFDX ,W ) 
IF(K0UNT)5000,115,114 

114  Y(1)=XX(KCUNT ) 
IPCSIT(1)=3 

CALL  NAP1 (X,S,Y,IPOSIT, 1,C) 
CALL  FUNCS(C,8,F,8) 
FF(KCUNT,1 )=F(6) 
FF(KCUNT,2)=F(5) 
GO  TC  113 

115  CC2DX3=CFCX( 1 ) 
IF(V21.LE.0.)CC2DX3=C. 
D^1DX3=DFCX«2 ) 

CV2DX 3= CC2CX 3/(1078. -Q.6*X(2>) 
DVODX3=DQ20X3/I1078.-0.6*X(4) ) 

0ERY(3)=-LAMCA(2)»( ( X ( 2  )-HlV  )  *DV2DX3  +  DQ 10X3 ) /X ( 1 ) +LAMDA ( 3 ) * 
1DV0DX3  +  LACDA(4)*( FEEC*(HF-X(4 )  J+Q2  +  V0  2M X(4)-H2V))/(X(3)+X( 

3)  )- 
2LAMDA(4)MCU2CX3-HX(4  )-H2V  )* C V0DX3 ) / X ( 3 )  +  (LAMDA(  1  )*0V2D 

X3) 
IF(X(3) .LT.l. )DERY(3 )=0. 
GC    TC    200 

ALL  VARIABLES  AT  STEACY  STATE 


179 


125 


5C0O 
5010 

200 


DERY ( 1 ) =0. 

0ERY(2)=0. 
DERY(3)  =  '3. 
DERY(A)=0. 

C" D  TC  200 

WRITE (6,5C10)K0UNT 

FORMATtSX^KOUNT*' , 14) 

RETURN 

CONTINUE 

RETURN 

ERROR  EXIT.   I  OUT  OF  RANGE.   STOP  EXECUTION 


50  WRITE(6f2C50)  (  I(J)t J  =  lf A) 

STCP 
2050  F0RMAT(5X,«LAMD0T.  I  OUT  OF  RANGE%AI3) 

END 
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C  CONVAR 

C      PURPCSE 

C         DETERMINES  THE  OPTIMAL  CONTROL  LAW  BY  MINIMIZATION  OF 

C         THE  HAMILTONIAN 

C      USAGE 

C         CALL  CCNVAR(TOUT,LAMDA) 

C      PARAMETERS 

C         TCUT    -  TIME,  MINUTES,  AT  WHICH  CONTROLS  ARE  TO  BE 

C  FCUND,  INPUT 

C         LAMDA   -  INPUT  VECTOR  CONTAINING  ADJOINT  VARIABLES  AT 

C  TIME  TOUT 

C      SUBPROGRAMS  REQUIRED 

C         SUBROUTINE  SEARCH 

C         SUBROUTINE  MAPI 

C         SUeRCUTINE  FUNCS 

C         THE  FOLLOWING  SUBPROGRAMS  ARE  USED  FOR  THE  MI NI MI ZAR I  ON. 

C         THEY  ARE  NCT  LISTEC  SEPARATELY  IN  APPENDIX  B 

C         SUBROUTINE  PARTI 

C         SUERCUTINE  PART2 

C         SUBROUTINE  MAPCWN 

C         SUBROUTINE  MAPUP 

C         SUBROUTINE  VA04A 

C 

SUBROUTINE  CONVAR ( TOUT ,L AMDA ) 

REAL  LAMCA(A) 

CIMENSICN  UK),E(4),VU),F(8),Dm,X(8),W(60),IPOSIT(4), 
1ZU),XL0U),XUP(4),IL0(4),IUPU),ICLS(4),DX(4) 

C0MM0N/DELT/DELX1,DELX2,CELX3,.DELX4,DELU1,DELU2,DELU3,DELU4 

COMMCN/CCNTRL/UK  120  ) ,  U2  (  120  )  ,.U3  (  120  )  ,  U4  (  120  ) 

C0MMCN/CATA1/NDATA1,TEMP(22),THALP(22),THALPV(22) , 
1ALAMCA(22),PRESS(22) ,RH0V(22 ) 

C  OMMCK/ STATES/ XI  {  120),  X2(  120)  ,-X  3(120  ),X4(  120) 

COMMCN/STEACY/XHAT( A ),XINIT( A ) , T FrHF , NORDER , U1MAX , U1MI N , 
1U2MAX,U2MIN,U3MAX, U3M  IN, U4MAX ,U4M IN 

COMMCN/I  \CEX/  I  (4) 

COMMON/ CNEW/U1NEW( 120),U2NEW( 120),U3NEW( 120) ,  U4NEWU20) 

COMMCN/TI^ER/TIME(120),ITIME,  I  MAX  ,.  I  COUNT 

CAT  A  ZERC,FRChG,?ISCALE,SC,MAXIT,.  IP,  ICON/  1  .E-  5  ,  0.  5  , 1  .  E  +  70  , 
11.  E-3, 103,1,1/ 
C 

C         INTERPCLATE  FOR  VALUES  OF  THE  STATE  AND  ADJOINT 
C         VARIABLES  AT  TIME  TOUT 
C 

CALL  SEARCH(TIME,X1,T0UT, I MAX , NORDER , X ( 1) ) 

CALL  SEARCH (T  IME,X2,T0UT,  I  MAX ^NORDER , X ( 2 )  ) 

CALL  SEARCh(TIME,X3,T0UT,  I  MAX  ,  NORDER , X ( 3 )  ) 

CALL  SEARCH (TIME,X A, TOUT, I  MAX .NORDER , X ( 4 )  ) 

CALL  SEARCH (TIME,U2, TOUT, I MAX ,NORDER , X ( 6 ) ) 

TT=TCUT 
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C 
C 

c 
cc 


IF(X( 1J.LE.0. )X( l)=O.0Ol 

IF(X(3).LE.O.  )X{ 3)=0.07 

DO  1  J=5,8 

IFCXl J).LT.O.  )X( J  )=0. 

DC  3  J=l,4 

I ( J)=0 

IF(X( J) .LT.XHATt J  )  )GG  TO  3 

I { J)=l 

CGNTINUE 

THE  FGLLOWING  VALUES  OF  THE  CONTROL  VARIABLES  WERE 
SET  BECAUSE  OF  THE  FEEDBACK  LAW 

X(8)=U4KAX 
X(5)=U1MAX 

IF ( I (3).EC.1)X(5)=U1MIN 
X(7)=U3MAX 

I  F  (  I  (1).EG.0)X(7)=U3MIN 
FEED=X(5) 
00  2  J  =  l,4 
ILC( J)=l 
IUP( J)=l 

IGOTG=I( 1)+2»I(2)+4*I(3)+8*I(  4)  +  l 

GO  TC  (5!?0, 5, 200, 5, 2C%.5, 200, 5, 200,200, 200, 5, 200, 200, 200, 
15)  ,IGCTC 

THIS  SECTION  OF  THE  SUBPROGRAM  IS  EXECUTED  IN  ALL 
STAGES  EXCEPT  WHEN  NO  STATE  VARIABLE  IS  AT  STEADY  STATE 
Tl  T  TF 


PUT  CONTROL  VARIABLES  THAT  HAVE  TO  BE  SEARCHED  OVER 
IN  VECTOR  Z 


Z(1)=X(7) 
Z(2)=X(8) 
Z(3)=X(5) 
N=3 

SPECIFY  LOWER  (XLO) 
CCNTRCL  VARIABLES 


AND  UPPER  (XUP)  BOUNDS  ON  THE 


XL0(l)=U3f IN 
XL0(2)=U4MIN 
XL0(3)=U1MIN 
XUP( 3)=U1MAX 
XUP( 1 )=U3MAX 
XUP(2)=U4MAX 


SET  UP  PARAMETERS  FOR  MINIMIZATION  ROUTINE  VA04A 
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DC  7  J=1,N 
E( J)=SC*XUPl J  ) 
RANGE=XUP( J)-XLC( J  ) 
ESC=RANGE*FRCFG/E(J) 
IF(ESCALE.GT.ESC)ESCALE=ESC 
7  CONTINUE 
JTEST=1 
10  IT=3 

SUERCUTINE  PARTI  PARTITIONS  THE  VARIABLES  IN  Z  INTO 
THCSE  WHICH  CANNOT  BE  SEARCHEC  OVER  AND  INTO  THOSE  WHICH 
CAN  BE  SEARCHEC  OVER.   IF  A  CONSTRAINT  IS  ENCOUNTERED 
IT  NUMERICALLY  EVALUATES  THE  KUHN-TUCKER  MULTIPLIER 
ASSOCIATED  WITH  ThE  PARTICULAR  CONTROL  VARIABLE  TO 
DETERMINE  WHETHER  THE  CONSTRAINT  IS  TO  BE  HELD  OR 
RELEASEC.   IT  ALSC  DETERMINES  WHEN  THE  SEARCH  IS 
COMPLETE  AFTER  ACCOUNTING  FOR  THE  CONSTRAINTS 


CALL  PARTI  (Z,E,XLC,XUP,ILC,  I  UP,- 1  CL  S,  N  ,  H, 
IF( JTEST-A)80,20,90 


JTEST,DX) 


SUBROUTINE  MAPDWN  MAPS  THE  UNCONSTRAINED  VARIABLES  IN  Z 
AS  DETERMINED  IN  PARTI  INTO  A  VECTOR  U 

20  CALL  MAPCWN(Z,E, ICLS  ,N,U,V,MU) 
IFCMU.EQ.01G0  TO  10 
GC  TC  (30,40,50,60),  IT 

SUBROUTINE  VA04A  IS  THE  CONJUGATE  GRADIENT  SEARCH 
PROGRAM  WHICH  SEARCHES  OVER  THE  VARIABLES  IN  U  TO 
MINIMIZE  THE  HAMILTONIAN  H 

30  CALL  VAG4A1 

GC  TC  70 
40  CALL  VA04A2 

GC  TC  70 
50  CALL  VA04A(U,V,MU,H, ESCALE, IP, ICON,MAXIT, IT,W) 

GC  TC  70 
60  CALL  VA04A4 

SUERCUTINE  MAPUP  RE-MAPS  THE  VARIABLES  IN  VECTOR  U  INTO 
VECTC3  Z 

70  CALL  MAPUPIZ, ICLS  ,N,U,MU) 

SUBROUTINE  PART2  CHECKS  WHETHER  THE  SEARCH  CAUSED  ANY 
CF  THE  VARIABLES  TO  VIOLATE  CONSTRAINTS 

CALL  FART2(Z,XL0,XUP,  ILO,  IUP,  ICLS  ,N,IT) 
80  CONTINUE 
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THIS  SECTION  EVALUATES  THE  HAMILTONIAN  H. 


IPOS 
IPOS 
IPOS 
CALL 
CALL 
00  8 
'  X(J) 
H12= 
V21  = 
V02  = 
H12  = 
G1  =  F 
Q2  =  F 
HIV= 
H2V  = 
H=LA 
1»(X( 
2X(6) 
IF(  J 
GO  T 


IT  (  1 

IT(2 
IT(3 
PAP 
FUN 
7  J  = 
=  C(  J 
F(l) 
F(2) 
F(3) 
F  K) 
(5) 
(6) 
F(7) 
F(8) 
NCA( 
5)»< 
-V02 
TEST 
0  10 


)=7 

)=8 

)=5 

1  (X,  8f.Zf  IPOSITtNiD) 

CS(C,B,F,8) 

1,8 

) 


2)*(W12»(H12-X(2) )+V21»(X(2)-HlV)+Ql)/X(l)+LAMDA(4) 
HF-X14) )+Q2+VC2*(X(4)-H2V) ) /X ( 3 ) +LAMDA ( 3 ) * ( X { 5 ) - 
) 

.EQ.4.AND.IT.NE.5)G0  TO  20 


CONVERGENCE  IS  OBTAINED.   THE  OPTIMAL  VALUES  OF  THE 
CONTROL  VARIABLES  ARE  STORED  IN  ARRAYS  IN  COMMON  BLOCK 

CNEW 

90  U3NEW(ITIPE)=Z(1) 
U4NEH(ITIME)=Z(2) 
UlNEW(ITiyc)=Z(3) 
U2NEW( ITINE)=X(6) 
RETURN 

THIS  SECTION  OF  THE  SUBPROGRAM  IS  EXECUTED  MEN  NO 
STATE  VARIABLES  ARE  AT  STEADY  STATE.   T  Tl 

500  Z(1)=X(6) 
Z(2)=X(8) 
N=2 
XL0(1)=U2MIN 

XL0(2)=UAVIN 

XUP( 1)=U2PAX 

XUP(2)=U4^AX 

DO  505  J  =  1,N 

EU)  =  SC*XUP(J) 

RANGE  =  XUP( J)-XLO(  J) 

ESC=RANGE»FRCHG/E(J ) 

IF(ESCALE.GT.ESC)ESCALE=ESC 
505  CONTIKUl 

JTEST=1 
510  IT=3 
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CALL  PARTKZ,  E,XLO,XLP,  ILC,  IIP,  ICLS, N,H,   JTEST.DX) 

511  CONTINUE 

IF (JTEST-4)  580,520, 600 
520  CALL  PAPCWMZ,.E,  ICLS  ,N,U,V,MU) 

IFIPU.EC.OGO  TO  510 

GO  TC  (530,540,550,560),  IT 
530  CALL  VA04A1 

GO  TC  570 
540  CALL  VA04A2 

GO  TC  570 
550  CALL  VAG4MU,V,MU,H,  ESCALE,  IP,  ICON,MAXIT,  IT,W) 

GO  TC  570 
560  CALL  VA04A4 

570  CALL  yAPUP(Z,ICLS  ,N,U,MU) 

CALL  PART2(Z,XL0,  XUP,  ILO,  IUP,  ICLS  ,.N,IT) 

571  CONTINUE 
580  IPCSIT(1)=6 

IP0SIT(2)=8 

CALL  N"AP1  (X,8,Z,  IPCSIT,N,D) 

CALL  FUNCS(C,8,F,8) 

OC  590  J=l,8 

590  XI  J)  =  C( J) 

591  CONTINUE 
W12=F(1) 
V21=F(2) 
V02=F(3> 
H12=F(4) 
Q1=F(5) 
C2=F(6) 
H1V=F(7) 
H2V=F(8) 

H=LAVCA(1  )*(W12-V21-X(7)  )      +L AMDA ( 2 ) * ( W12* ( H12-X ( 2 ) )  +  V21 

*(X(2) 

1-H1V)+Q1)/X(1)+LAMCA(3)*(X(5)-XI 6  )-V0  2 ) +  LAMDA ( 4 ) * ( X ( 5 ) * ( HF- 

X(4)  ) 

2+Q2+VC2«(X{4)-H2V))/X(3) 

592  CCNTINUE 

IF( JTEST.EC.4.AND.IT.NE.5)  GO  TO  520 
GO  TG  51? 
600  U2NEW(ITIPE)=Z(1) 
U3NEW( ITINE )=C. 
U4NE*(ITI^E)=Z(2) 
U1NEW(ITIME)  =  X.I5) 
-tETURN 

ERROR  EXIT.   I  OUT  OF  RANGE.   STOP  EXECUTION 


200  WRITE(6,2050)  (  K J  ), J  =  l,4) 
2050  F0RVAT(2X,«  IN  CONVAR.   I  =  •  , 4 (  2X,  I  2  )  ) 
STOP 
END 
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C  INPUT 

c***« *****»«*«*******»»««*«*************** *********************** 


PURPOSE 

TC  INPUT  DATA  FOR  COMMON  BLOCKS  DATA1,  DATA2,  STEADY 

AND  CELT 
USAGE 

CALL  INPUT 
SUBPROGRAMS  REQUIRED 

SUBRCUTINE  SEARCH 

SUBROUTINE  INPUT 

C0MMCN/CATA1/NCATA1,TEMP( 2  2),THALP( 22 ) , THALP V ( 22 ) , 
1ALAMCA(22),PRESS( 22),RHOV(22) 

COMMON/ ST EACY/XHAT( 4  ),XINIT( 4  )  , TF , HF , NOROER , U1MA X , Ul Ml N, 
1U2MAX,U2MIN,U3MAX,U3MIN,U4MAX,UAMIN 

C0MMCN/DELT/DELX1,DELX2,0ELX3,DELX4,DELU1,DELU2,DELU3,DELU4 

COMMCN/DATA2/NDATA2»T( 8) ,  RHO( 8 ) ,CP ( 8 ) , VI SC I  8 ) ,THCOND(8)  t 
1PRANTL(8),BETA 

READ  NUMBER  OF  POINTS  IN  ARRAYS  IN  DATAU  AND  DATA2, 
CRCER  OF  INTERPOLATING  POLYNOMIAL,  FEED  TEMPERATURE 
ANC  EXPANSION  COEFFICIENT  OF  WATER 

READ(5,10CO)NCATA1,NCATA2,NORCER,TF,.BETA 

READ  PROPERTIES  OF  WATER  AND  STEAM  VS  TEMPERATURE,  I.E. 
ENTHALPY, VAPOR  ENTHALPY,  LATENT  HEAT,.  SATURATION 
PRESSURE  AND  VAPOR  DENSITY 

READ (5, 1010) (TEMPI  I)  ,THALP( I  ), THALP V( I ),ALAMDA<I ) .PRESS (I ) , 
1RH0V( I )  ,  I  =  1,NCATA1) 

READ  PROPERTIES  OF  WATER  VS  TEMPERATURE,  I.E.  DENSITY, 
SPECIFIC  HEAT,  VISCOSITY,  THERMAL  CONDUCTIVITY  AND 

PRANDTL  NUMBER 

READ(5,  1010  MTU  )  ,RHC(  I  ),CP{  I  ),VISC(  I),THCOND(I)  ,PRANTL(I)  , 
1I=1,NCATA2) 

FIND  FEED  ENTHALPY  CORRESPONDING  TO  FEED  TEMPERATURE 
CALL  SEARCH(TEMP,THALP,TF,NDATA1,N0RDER,HF) 

READ  STEADY  STATE  VALUES 
READ(5,10  30)XHAT( 1 ) , T1SS , XHAT ( 3 ) , T2SS , X3MAX 

DETERMINE  STEADY  STATE  ENTHALPIES  AND  TEMPERATURES 
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CALL  SEARCH < T EMP , THAL P , T 1 SS , N DATA  1 , NORDER , XHAT ( 2 )  ) 

CALL  SEARCMTEMP,THALP,T2SS,NDATA1,N0RDER,XHAT(4)  ) 
C 

C         READ  INITIAL  VALUES  OF  STATE  VARIABLES 
C 

READ(5,103C) (XINITl  I ), I  =  U4) 
C 

C         READ  MAXIMUM  AND  MINIMUM  VALUES  OF  CONTROL  VARIABLES 
C 

READ(5,1040)U1MAX,U1MIN,U2MAX,U2MIN,U3MAX,U3MIN,U4MAX,U4MIN 

C 

C         READ  INITIAL  PERTURBATIONS  IN  STATE  AND  CONTROLS  FOR 

C         USE  IN  SUBROUTINE  PCERIV 

C 

REACC5,1040)DELX1,CELX2, CELX  3  ,DELX4,  DELU  1  ,.DELU2,  DELU3  ,  DE  LU4 
C 

C  CUTPUT    DATA    READ     IN 

C 

WRITE(6,2000) 

WRITE(6,2010)(TEMP(I),THALP(I),THALPV(I) .ALAMDAd )  , 
1  PRESS ( I ) , RHOV( I),I=1,NCATA1) 
WRITt (6,2020) 

WRITE(6,20  30)(T(I),RI-C(I),CP(I  ),THCOND(  I  )  ,  V  I  SC  ( I  )  , 
1PRANTH  I  )  ,  I  =  1,NDATA2  ) 
WRITE(6,2040)TF, HF,BETA 

WRITE (6,2050) ( I, X INIT( I ) ,XHAT( I), 1=1,4) 
WRITE(6,2080) 

WRITE(6,2C9U)U1MAX,U1MIN,U2MAX,U2MIN,U3MAX,U3MIN,U4MAX, 
LU4MIN 

WRITE(6,3040)CELXl,CELX2,DELX3,DELX4,DELUl,DELU2,DELU3, 
1DELU4 
RETURN 
1000    FORMAT(3I2,3F10.0) 
1C1C    F0RMAT(6F10.0) 
1020     F0RMATU2) 
1C30    FCRMAT(5F10.0  ) 
1040    F0RMAT(8F10.     ) 

1C60    FURMAT(  IX  ,  I  2  ,  2X  ,  F6  .  3  ,  F  10  .  7,.F  10  .  6  ,  F10.  7,  F  10.6  ,  3F8  .  5,  F  5  .  1 ) 
2000    FORMAT( lhl, 1CX,« PROPERTIES    OF    WATER    AND    STEAM •,//, 5X ,• TE MP • 
1,2X,'LIC    ENTHALPY' ,2X, «VAP    ENTHALPY •, 2X,  '  LATENT    HE A T •,  3X  , 
2«  PRESSURE*  ,2X,  «VAP    D  ENS  I  T  Y  •  ,  /  ,.4X  ,  •  (  DEG    F  )  •  ,  3X  ,  •  (  B  TU/LB  )  •  ,  5X 
B.MBTU/LBPfTXtM  8TU/LB)  •  ,  5X  ,  •  (  PS  I  )  •  ,  4X ,  •  (LBS/CU    FT)'  J 
2010    F0RMAT(4X,F5.1,4X,F6.2,7X,F7.2,.8X,F7.2,5X,F7.4,4X,F7.5) 
2  020    FCRMAT(//5X,'TEMP«,3X,'DENSITY,,4X,«SP.HT.  •  ,  2  X  ,  •  TH.COND.  '  ,• 
13X,'VISCCSITY' ,2X, 'PRANDTL    NO  •  , /  ,  4X ,  • ( DEG    F)',1X, 
2,(LBS/CUFT),,1X,MB/LB.F)«,1X,MBF/HSQFF)«,1X,  MLBS/FTHR)  •  ) 
2030    F0RMAT(4X,F5.  1,4X,F5.2,5X,F6.4,.5X,F5.3,4X,F6.4,5X,F5.3) 
2040    FORMAT! lhl,5X, 'STEADY    STATE    AND     INITIAL    VALUE S  ',/, 5X , 

l'FEED    TEKPERATURESFS.lt/iSX,  'FEED    ENTHALPY  =  •  ,  F  7.  2  ,  /  ,  5X  , 

2'BETA',12X,'=',F9.7,//8X,'XINIT',6X,'XHAT') 
2 0  50    FORMAT(5X,Il,2X,F5.2,5X,F6.2) 
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2C80  F0RMAT(1H1,5X, 'MAX  AND  MIN  VALUES  OF  CONTROL  VARIABLES',//, 

117X,'FAX  VALUE', 5X, 'fIN  VALUE') 
2C90  FORMAT(18X,F7.3,7X,F7.3) 
3040  FORMAT(1H1,10X,' INITIAL  INCREMENT  IN  INDEPENDENT  VARIABLES' 

1,//,5X,'CELX1  =  ',F5.2,/,.5X,  • DELX2= • ,F 5. 2 , / , 5X ,  •  DELX3=  •  ,  F5  .  2  , 

2/,5X,'DELX4=',F5.2,/,5X,'DELUl=',F5.2,/,5X,'DELU2=',F5.2,/, 

3  5X,' CELU3=',F5.2,/,5X,'DELU4=',F5.2) 
END 


C  SEARCH 


PURPCSE 

TO  INTERPOLATE  USING  A  POLYNOMIAL  OF  DESIRED  ORDER 

USAGE 

CALL  SEARCH! X,Y,XG I VEN , NDATA , NORDER , YFOUND > 

DESCRIPTION  OF  PARAMETERS 

X      -INPUT  VECTOR  OF  INDEPENDENT  VARIABLE 
Y      -INPUT  VECTOR  OF  DEPENDENT  VARIABLE 
XGIVEN-INPUT  POINT  AT  WHICH  VALUE  OF  Y  IS  REQUIRED 
NCATA  -INPUT  NUMBER  OF  CATA  POINTS 
NORCER-INPUT  ORCER  OF  INTERPOLATING  POLYNOMIAL 
YFCUND-OUTPUT  VALUE  OF  CEPENDENT  VARIABLE 

REMARKS 

MAX  CIMENSION  OF  INPUT  VECTOR  IS  20  AND  MAX  ORDER  OF 
INTERPCLATING  PCLYNOMIAL  IS  4  (IF  GREATER , CHANGE 
CIMENSION  STATEMENTS  ACCORDINGLY 

SUBPRCGRAMS  REQUIRED 

FUN CT  I  CM  AITKEN(X,Y,NDATA,XGIVEN) 

CcTHCC 

TAKES  NORCER+1  POINTS  AND  INTERPOLATES  USING  AITKEN- 
LAGRANGE  INTERPOLATION 


SUBROUTINE  SE ARCH ( X, Y, XG IV  EN , NDATA, NORDER , YFOUND ) 
DIMENSION  X(20),Y(20),A{5),B(5) 

SEARCH  FOR  NORDER  POINTS 

YFCUNC=0. 

DC  30  M=l, NCATA 

IF(XGIVEN-X(M)  >10,20,30 
10  JJ=N 

GO  TC  40 
20     YFCU.\C  =  Y(M) 

RETURN 
30  CONTINUE 
40  LL=JJ-NCRCER/2+l 

IF(LL.LE.0)LL=1 

LU=LL+NCRCER-1 

IFILU.LT. NCATA)  GO  TC  50 

LU=NCATA 

LL=LU-NCRCER+1 
50  DO  60  M=LLtLU 

K=M-LL+1 

A(K)=X(M) 

BKK)=Y(M) 
60  CONTINUE 

IF(YFCUNC.EC.0)YFOUNC  =  AITKEN( A, B , NORDER, XGI VEN  ) 

RETURN 
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END 
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£•«•««««»»•«»«•«*«««*•*•«*#»•*«•*##••*•**»«••*»»****«*•»**•*•*»*» 

C  AITKEN 

c    purpc.se 

C  INTERPOLATES  USING  POLYNOMIAL  OF  DESIRED  ORDER 

C      USAGE 

C  FUNCTION  AITKEN(XtY,NCATA,XGIVEN) 

C      DESCRIPTION  OF  PARAMETERS 

C  X      -I!\PUT  VECTOR  OF  INDEPENDENT  VARIABLE 

C  Y      -INPUT  VECTOR  OF  CEPENDENT  VARIABLE 

C  NCATA  -NUMBER  OF  DATA  POINTS 

C  XGIVEN-INPUT  POINT  AT  WHICH  VALUE  OF  Y  IS  REQUIRED 

C      REMARKS 

C  XGIVEN  MUST  LIE  IN  THE  RANGE  OF  THE  TABLE.  IF  NDATA  IS 

C  GREATER  THAN  5,  ADJUST  CIMENSION  STATEMENT 

C      METHCC 

C  AITKEN-LAGRANGE  INTERPOLATION 

C 

FUNCTION  A  I TK EN (X,Y, NCATA, XGIVEN  ) 

DIMENSION  X(5),Y(5),Z<5,5) 

DC  10  1=1, NCATA 
10  Z(I,1)=Y(  I) 

NN=NCATA-1 

DC  20  L=1,NN 

LL=L+1 

DC  30  K=LL, NDATA 
30  UK, L+l)=( (X{K)-XGIVEN)»Z(L,L)-(X(L)-XGIVEN)*Z(K,L) ) /(X(K)- 

X(L)  ) 
20  CONTINUE 

AITKEN=Z(NDATA, NDATA) 

RETURN 

END 
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C#**«# ***»»*«*«*»*»««************* ******************************* 

C  FUNCS 

c**»**«»»«******************************************************* 

C  PURPOSE 

C  CALCULATES  VALUES  OF  ALGEBRAIC  VARIABLES  FROM 

C  VALUES  OF  CONTROL  ANC  DIFFERENTIATED  VARIABLES. 

C  USAGE 

C  CALL  FUNCS(X,N,.F,M) 

C  DESCRIPTICN  OF  PARAMETERSO 

C  X  -  N  DIMENSIONAL  VECTOR  OF  DIFFERENTIATED  AND 

C  CONTROL  VARIABLES  AS  DESCRIBED  IN  SUBROUTINE  XDOT 

C  F  -  M  CIMENSIONAL  VECTOR  CONTAINING  VALUES  OF 

C  ALGEBRAIC  VAR IABLESt WHERE, 

C  F(l)  -  FEED  TO  FIRST  EFFECT,  W12  (LBS/MIN) 

C  F(2)  -  VAPOR  FROM  FIRST  EFFECT,  V21  (LBS/MIN) 

C  F(3)  -  VAPOR  FROM  SECOND  EFFECT,  V02  (LBS/MIN) 

C  F(4)  -  ENTHALPY  OF  FEED  TO  FIRST  EFFECT,  H12,. 

C  BTU/LB 

C  F(5)  -  HEAT  TRANSFER  IN  FIRST  EFFECT,  Ql, 

C  OTU/MIN 

C  F(6)  -  FFAT  TRANSFER  IN  SECOND  EFFECT,  Q2 

C  BTU/MIN 

C  F(7)  -  VAPOR  ENTHALPY  IN  FIRST  EFFECT,  HIV, 

C  BTU/LB 

C  F(8)  -  VAPOR  ENTHALPY  IN  SECOND  EFFECT,  H2V, 

C  BTU/LB 

C  SUBPROGRAMS  REGUIRED 

C  SUBROUTINE  H^AT 

C  SUBROUTINE  SEARCH 

C 

C  IN  THE  RANGE  T  =  212  DEC  F  TO  220  DEG  F  (I.E.  ENTHALPY 

C  180.07  BTU/LB  TO  188.13  BTU/LB) 

C  RELATIONSHIP  BETWEEN  VAPOR  ENTHALPY  AND  TEMPERATUREO 

C  HV  =  1G66.0  +  O.A*T 

C  RELATIONSHIP  BETWEEN  TEMPERATURE  AND  LIQUID  ENTHALPYO 

C  T  =  HL  +  32.0 

C  WHERE, 

C  T  =  TEMPERATURE  (CEG  F) 

C  HL  =  LIQUIC  ENTHALPY  (BTU/LB) 

C  whERE  hV   =  VAPOR  ENTHALPY  (BTU/LB) 

C  RHO  =  VAPOR  DENSITY  (LBS/CU.FT) 

C 

DIMENSION  X(l  )  ,F(  1) 

C0MMCN/CATA1/NDATA1,TEMP( 22 ) , THALP ( 22 ) , THALP V( 22 ) , 
1ALAMDA(22  ), PRESS ( 22) ,RHOV( 22) 

COMMON/ST  EADY/XHAT(  4  ),X  I  NITU),TF,HF,N0RDER,U1MAX,U1MIN, 

1U2MAX,U2MIN,U3MAX,U3MIN,U4MAX,U4MIN 
CONMCN/INCEX/ I (4) 
C 
C         FIND  TEMPERATURES  AND  VAPOR  ENTHALPIES  IN  EFFECTS 

C 
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5  Tl=X(2)+32.0 

T2=X(4)+32.0 

H1V=1066.+0.4*T1 

H2V=1066.+0.4*T2 

I  F  (  I  (2).EUi)G0  TO  10 
C 

C         X2.NE.XFAT(2),  FIRST  EFFECT  IS  BEING  HEATED. 
C 

W2=0. 

V21  =  . 

VC2=0. 

M2=X(6)+X(7) 

IF(W12.NE.0.)G0  TO  1 

T12=T1 

GC  TC  2 

1  H12=(X(7)«X(2 )+X(6)*X(4)  J/W12 
T12=H12+32. 

2  CALL  HEAT(CltX(l),W12iV21.X(8)i.TlfT12»l,i) 
GC  TC  70 

C 

C         X2.EU.XHAT(2) ,  FIRST  EFFECT  IS  BOILING. 

C 

10  CONTINUE 

CALL  FEAT(C2,X(3),0.,0.,T1,T2,TF,3,2) 

CALL  SEARCH { T EMP, ALAMnAt Tlf NDATA It NORDERtAL AMI) 

V21=G2/ALAM1 

IF(I (1).EC.1)X(6)=V21 

W12=X(6)+X(7) 

IF(W12.NE.C.)G0  TO  190 

T12=T1 

GC  TC  20  0 
190  T12=(X(7)*TH-X(6)*T2)/W12 

H12=T12-32. 
200  CALL  FEAT(G1,X(1 ) ,W12,V21,X< 8)„Tl,T12, 2, 1) 

I F { I (4).EC.0)G0  TO  70 
C 

C         X4  =  XHATU),  SECONC  EFFECT  IS  ALSO  BOILING. 
C 

CALL  SEARCH (TEMP, AL ANTA , T 2, NDATA 1 , NORDER , ALAM2 > 

V02=(C2-X(5)*(X(4)-HF))/ALAM2 

GO  TO  80 
70  VC2-  . 
30  F(1)=U12 

F (2)=V21 

F(3)=V02 

F(4)=H12 

F(5)-tl 

F(6)=G2 

F17)=H1V 

F(8)=F2V 

kETURN 
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C*** 

c 
c»  *« 

c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 


♦♦•••♦••••••••♦•A******************************************** 

HEAT 


THE    HEAT    TRANSFER    RATE 


PURPCSE 

EVALUATES 
USAGE 

CALL  HEAT(C«H,W,V,TO,TI,TIN, IPHASE.IEFECT) 

PARAMETERS 
Q 


H 

W 

V 

TO 

TI 

TIN 

IPHASE  - 


IEFECT  - 

SUBPROGRAMS 
FUNCTION 
FUNCTION 
FUNCTION 
FUNCTICN 
FUNCTION 


HE 
HO 
LI 
VA 
OU 
IN 
IN 
PA 


PA 

Rl 
FC 
FC 
FC 
FC 
HE 


SUBROUTINE 


AT  TRAN 
LO-UP, 
GUIC  FL 
POR  FLO 
TSICE  ( 
SIDE  (L 
LET  TEM 
RAMETER 
It  INCI 

2,  INCI 

3,  INCI 
RAMETER 
ANSFER 
GUI  RED 
INB 

IN 

0 

IF 

IGHT 

ALPHA 


BTU/MIN,  OUTPUT 


SFER  RATE, 

LBS,  INPUT 

OW  RATE,  LBS/HR,  INPUT 

W  RATE,  LBS/HR,  INPUT 

STEAM  SICE)  TEMPERATURE,  DEG  F,  INPUT 

IQUID  SIDE)  TEMPERATURE,  DEG  F,  INPUT 

PERATURE  OF  LIQUID,  DEG  F,  INPUT 

DESCRIPTIVE  OF  TYPE  OF  BOILING 
CATES  SINGLE  PHASE  FLOW 
CATES  TWO  PHASE  FLOW 
CATES  NATURAL  CONVECTION 

INDICATING  EFFECT  FOR  WHICH  HEAT 
RATE  IS  TO  SE    CALCULATED,  INPUT 


SUBROUTINE  HE  AT  (  Q,  H,  W,  V,  TO,  T  I  ,.T  IN,  IPHASE,  IEFECT) 

REAL  LIS, LI, LIB 

DIMENSION    THETA{ 3) 

COMMON/ C  AT  A 1/ N  DAT  A  1,  TEMP  (22)  ,  THALP  (  22  )  ,.THALP  V  (  22  )  , 
1ALAMDA(22) , PRESS ( 22 ) , RHOV ( 22 ) 

COMMCN/DATA2/NCATA2,T(8),RHO( 8),CP(8),VISC{8) ,THC0ND(8) , 
1PRANTH8  )  ,BETA 

COMMCN/INCEX/ I (4) 

DATA  NCRCER/4/ 

IFUC.LE.TI  )GO  TO  25 

LI  IS  THE  MAXIMUM  LENGTH  OF  THE  TUBES  IN  FEET. 
A  IS  THE  MAXIMUM  HEAT  TRANSFER  AREA  IN  SQ.FT. 

Ll=9.5 

IF(IEFECT.EC2)L1  =  1.917 

-=7.46128 

IFCW.EQ.O. )TIN=TI 

TBULK=(TII\  +  TI  )*0.5 

IFtV.LE.}. )IPFASE=1 

IF(  IEFECT. EQ. 2) A= 9. 08116 
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10 


11 
12 


THE  CUTER  WALL  TEMPERATURE,  TW2  AND  INNER  WALL 
TEMPERATURE  TWl  ARE  GUESSED  INITIALLY  AND  THEN  ITERATED 
ON  BY  EVALUATION  CF  THE  FILM  COEFFICIENTS  AND  THEN 
RE-EVALUATING  THE  WALL  TEMPERATURES 

TW2=G.5»(TC+TBULK) 
TW1=TW2 

ITIME=ITIME+1 

THE  PARAMETERS  THETA  ARE  CORRECTION  PARAMETERS  FOR  THE 


TWC  EFFECTS 

THETA1  =  0.1395 

=  0.0928 

THETA2  =  0.1359 

=  0.1763 


WHEN  FIRST  EFFECT  IS  NOT  BOILING 

OTHERWISE 

WHEN  SECOND  EFFECT  IS  NOT  BOILING 

OTHERWISE 


THETA(1)=0.1395 

THETA(2)=1. 

THETA(3)=1. 

I  F ( I (2) .EG.OJGO  TO  5 

THETA(1)=0.0928 

THETA(2)=0.1359 

THETA(3)=0.1763 

CONTINUE 

KCUNT=1 

TH1=TFETA(1) 

IFUEFECT.EQ.UGO    TO    11 

TH1=THETA(2) 

IF(I(4).EC.1)TH1=THETA(3) 

GC    TC    12 

IF(W.NE.O. )G0    TO    15 

IF ( I (4).NE.1.0R.IEFECT.NE.2)G0    TO    13 

HIS  IS  THE  INSIDE  FILM  COEFFICIENT 

HCS  IS  THE  OUTSICE  FILM  COEFFICIENT 

TWIN  AND  TW2N  ARE  THE  NEW  INNER  AND  OUTER  WALL 

TEMPERATURES 


HIS=FCING(TWlt.TI  ) 

GO  TC  17 
13  HIS=FCIN(TBULK,TW1»H, IEFECT) 

GC  TC  17 
15  HIS=FCIF(TEULK,W,TW1) 
17  H0S=FC0(TC,TW2, IEFECT) 

HW=849.62A 

TW1N=(TBULK*TH1  *H IS  * ( HOS  +  HW  )  +HW»HOS* TO ) / ( TH1  *HIS» 

1 (HOS+HW)+hOS»HW) 

TW2N=(H0S«T0  +  HW*TW1N  )/( HOS+HW) 

IF(ABS(TW1N-TW1).LT.1..AND.ABS(TW2N-TW2).LT.1. )G0  TO  20 

IF{KCUNT.GT.15)WRITE(6,21CO)HIS,HOS,TW1,TW2,TW1N,TW2N 

IF(IU) .NE.1.CR.IEFECT.NE.2.0R.KOUNT.GT.10)GO  TO  18 
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TWl=Twl+J.2*(TWlN 
TW2=TW2+0.2*( TW2N 
GC    TC    19 

18  TW1=TW1N 
TW2  =  TW2i\ 

19  KOUNT  =  KCUi\T  +  l 
IFIKCUNT.LT.20  )GO 
WRITE(6,2000)KOUN 
STOP 

20  IF(ABS(HIS).LT..1E-70)G0  TO  25 


l-TWl  ) 
I-TW2) 


TO  LO 

T,Tl.,TI,TIN,IPHASE,IEFECT 


CALCULATE  OVER 
SENSIBLE  HEATI 
LIS  IS  THE  LEM 
DELT1S  IS  THE 
SENSIBLE  HEATI 
CIS  IS  THE  HEA 
ZCNE 


ALL  HEAT  TRANSFER  COEFFICIENT  UlS  IN 

NG  ZONE 

GTH  OF  THE  SENSIBLE  HEATING  ZONE 

MEAN  TEMPERATURE  DIFFERENCE  IN  THE 

NG  ZONE 

T  TRANSFER  RATE  IN  THE  SENSIBLE  HEATING 


25 


U1S=(1./(1. 25357/ 

l1s=height( iphas5 
deltis=tc-(ti+;tin 
q1s=u1s*a  «l1s»de 
if( ipfase.eq.2)g0 

C  =  G1S 

RETURN 

Q=0- 

RETURN 


(TH1      *HIS)+1./HOS+0.001315+0.002) )/60. 
,Hf  TBULK,TW1,TIN,.W,IEFECT) 
)*0.5 
LT1S/L1 
TO  30 


THIS  PART  OF  T 
CALCULATION  OF 
EFFECT  WHEN  BO 
THE  EQUATIONS 


HE  SUBROUTINE  IS  USED  ONLY  FOR  THE 

THE  HEAT  TRANSFER  RATE  IN  THE  FIRST 
RING  OCCURS.   REFER  TO  APPENDIX  A  FOR 


30  TW2  =  u.5»(TC+TI  ) 
T  In  1  =  T  V.  2 

I TI ME= IT  INE  +  1 

31  CONTINUE 
TB  =  TI 
KCUNT=1 

32  IF(  (W-V) .GT.O.  )GQ  TO  35 
X12=l. 
HIB=22.15727*V«*0.8 

GC  TC  AO 
35  X12=V/W 

HIB=FCIF(TB,W  ,TW1  ) 
40  HCB=FC0(TC,TW2, IEFECT  ) 

HW=849.624 

TW1N  = (T8*THETA(1  ) 
1HIB*(H0B+HW)+H0B* 

TW2N  = (HCB»TC+HW#T 


*HIE«(HOB+HW )+Hfc«HOB*TO )/( THETA(l)* 

HU  ) 

WIN  )/(HOB  +  HW  ) 
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IF( ABS(TwlN-TWl) .LT. 1 .  .  AND.  ABS( TW2N-TW2 ) .LT.l.     )G0    TO    50 

TW1=TW1N 

TW2  =  TVv2N 

IF(TW1.LT.TI)TW1=TI 

IF(TW2.GT.TO)TW2=TO 

KCUNT=KCUNT+1 

IF(KCUNT.LT.10)G0  TO  32 

wRITE(6,2C00)K0UNT,TC,TI,TIN, I  PHASE, I EFECT 

STOP 
50  IF( (W-V) .LE.O. )G0  TO  65 

X12PR=0.4«X12 

CALL  SEARCH(T,RHO,TB    , N0ATA2, NORDER ,RHOL ) 

CALL  SEARCH(T, VISCtTB    , NDATA2, NORDER, VI SCL ) 

CALL  SEARCh(TEMP,RHOV,TI,NDATAl,NORDER,RHOG) 

VISCG=0.0195  3  29+5.10  7E-0  5*TI 

XTT=( (l.-X12)/X12 )**0.9*(RHOG/RHOL)**0.5*(VISCL/VISCG)**G.l 

XTTPR=( ( l.-X12PR)/X12PR)**O.9*(RHOG/RH0L)**0.5*( VISCL/VISCG 
1 >  **0- 1 

HITP  =  3.5*HIB/SGRTUTT) 

GT=4844.673*W 

CALL  ALPHA(XTT,GT,ALFA12) 

CALL    ALPI-A{XTTPR,GT,  ALFAPR) 

ALFA^=0.5*( ALFA12+ALFAPR) 

HINB=0.074*(TW1-TI )*«2.86 

HIB=ALFAM*HINE+HITP 
65  U1B=( l./( 1.25357/ (ThETA( 1 ) *H I B )+ 1 . /H0B+0.C01 315+0. 002 ) )/60. 

LlS=HEIGhT(  IPhASE,H,TBULK,TWl,TlN,W,.IEFECT) 

L1B=L1-L1S 

DELT1S=TC-(TI+TIN)*0.5 

Q1S=U1S*A  «L1S*DELT1S/L1 

DELT1E=TC-TI 

Q1B=U1B*A  *L1B*DELT1B/L1 

G=Q1S+Q1B 

RETURN 
2000  F0RMAT12X,' AFTER' , I4,2X, • ITERATIONS  IN  HEAT  THERE  WAS  NO  ', 
1' CONVERGENCE' ,/,2X, ' ARGUMENTS  WERE',/,5X, • TO = • ,F 10. 3 , 2X , 
2'TI=' ,F1J.3,2X,'TIN=',F10.3,  'PHASE=f , I4,2X,  • EFFECT= • , 14 ) 

END 


198 


C  HEIGHT 

C»**« •«««*«****««**«««**•*«*«********* *************************** 

C      PURPOSE 

C         DETERMINES  LIQUIC  HEIGHT  IN  EACH  EFFECT 

C      USAGE 

FUNCTION  HEIGHT  (  I  PHAS  E»  H,  TX  ,  T  1BW  ,  T  IN,  FLOW  ,  I  EFEC  T  ) 
C      PARAMETERS 

C         IPFASE   -  =1,  SINGLE  PHASE  FLOW 
C  =  2,  TWO  PHASE  FLOW 

C         H        -  HOLD  UP  OF  LIQUIC  (LBS) 
C         TX       -  MEAN  LIQUID  TEMPERATURE  (DEG  F) 
C         T1BW     -  INNER  WALL  TEMPERATURE.  LIQUID  SIDE.  (DEG  F) 
C         TIN      -  ENTERING  LIQUID  TEMPERATURE  (DEG  F) 
C         FLOW     -  FLOW  RATE  OF  ENTERING  LIQUID  (LBS/MIN) 
C         IEFECT   -  =1,  FIRST  EFFECT 
C  -  =  2,  SECCNC  EFFECT 

C      SUBPROGRAMS  RECUIREC 
C         SUBROUTINE  SEARCH 

C         SUBROUTINE  DET5  (FROM  THE  IBM  SSP ) 
C 

FUNCTION  HEIGHTl I PHA S E , H , TX, T 1BW , T IN , FLOW ,1 E FECT ) 

DINENSIGN  Z(20) 

COMMON/ C AT Al / N CAT  A  1,  TEMP  (22) ,THALP( 2 2 ) , THALP V ( 22 )  » 
1ALAMCM22),PRESS(22)  ,RHOV(22) 

COMMGN/CATA2/NCATA2,T(8),RHO(8),CP(8),VISC(8) ,THCOND(8) , 
1PRANTH8)  ,BETA 

CATA    NORDERM/ 

AL=2.3125 
C 

C         TUBE  CIAMETER  IN  BOTH  EFFECTS  =  0.0725  FT. 
C 

C=0.0725 
C 

C         NUMBER  OF  TUBES.  FIRST  EFFECT  =  3.   SECOND  EFFECT  =  15. 
C 

NT  =  3 

IF(  IEFECT. EC. 2)NT  =  15 

TBULK=TX 

CALL  SEARCH(T,RH0,TBULK,NDATA2,N0RDER,RH0F) 

IF(  IPHASE. EQ.2JG0  TO  20 

AT=3.14159*C»CM. 
C 

C         DEAD  HCLDUPS  (DO  NOT  CONTRIBUTE  TO  HEIGHT). 
C  FIRST  EFFECT  =  1.5315  LBS.   SECOND  EFFECT  =  16.336  LBS. 

C 

HCLC  =  H-1.5315 

IF(IEFECT.EQ.2)H0LD=H-16.3  36 

HEIGHT =HCLC/( FLOAT ( NT ) »RHOF* AT ) 

IFlHEIGHT.Lr. 0. )HEIGFT=0. 

IF(  IEFECT. EQ.2JG0  TO  10 
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LENGTH  CF  TUBES.  (ALSO  MAXIMUM  HEIGHT) 

FIRST  EFFECT  =9.5  FT.   SECOND  EFFECT  =  1.917  FT. 

IF(HEIGHT.GE.9.5)HEIGFT=9.5 
RETURN 

THIS  SECTION  OF  THE  SUBROUTINE  DETERMINES  THE  LENGTH  OF 
THE  SENSIBLE  HEATING  ZONE  IN  TWO-PHASE  FLOW.   REFER 
TO  APPENDIX  A  FOR  THE  EQUATIONS  AND  REFERENCES 

10  IF(HEIGHT.GE.1.917)HEIGHT=1.917 

RETURN 
20  HL=FCIF(TBULK,FLOW) 

CALL  SEARCH(T,CP,TBULK,NDATA2,N0RDER,CPB ) 

IF(FLCW.EC.O.  )G0  TO  30 

DTDL=(3.14159*D*FL0AT(NT)*HL*(T1BW-TBULK) ) / ( 60.*FLOW»CPB ) 
30  CPDL=-RH0F/144. 

CALL  CET5(10.»PRESS,Z,NDATAlt  IER ) 

CALL  SEARCH (TEMP,Z,T BULK, NDATA1 , NORDER, DPDT ) 

DTDP=1./DPDT 

IF(FLCW.EC.O. )G0  TO  40 

HEIGHT=AL«CTDP/     (-( CTCL/CPDL  MDTCP) 

GC  TC  50 
40  DELT=2.*(TX-TIN) 

HEIGHT=1.+(CELT/(DPDL*DTDP) ) 
50  IF(HEIGHT.GT.9.5 )HEIGHT=9.5 

RETURN 

END 
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FCINB 


PURPCSE 

CALC 

CCNV 
USAGE 

FUNCTICN  FCINB(TW,TI) 
PARANET 

TI 
SUBPRCG 
SUER 


ULATES  THE  INSIDE  FILM  COEFFICIENT  FOR  NATURAL 
ECTION  AND  POOL  BOILING  USING  THE  RHOSENOW  EQUATION 


ERS 

-  INNER  WALL  TEMPERATURE,  DEG  F,  INPUT 

-  INNER  BULK  TEMPERATURE*  DEG  F,  INPUT 
RAMS  REQUIRED 
CUTINE  SEARCH 


FUNCTION  FCINB(TW,TI  ) 

COMMON/CATA2/NDATA2,T<8) , RHO < 8 ) , CP < 8 )  ,  VI SC ( 8  )  , THCOND ( 8  )  , 
1PRANTH8)  ,8ETA 

COMMON/ DATA  1/ NDAT Al, TEMP ( 22) , THALP ( 2  2 ) , THALP V ( 22 ) , 
1ALAMDA(22),PRESS(22),RH0V(22) 

DATA  NCRCERM/ 

DELTAT  =  TVi-TI 

IFIDELTAT.LE.O. ) GO  TC  10 

TFILM=TW-0.75*(TW-TI  ) 

CALL  SEARCF(T,VISC,TI,NCATA2,NORCER,VISCF) 

CALL  SEARCF(TEMP,THALPV,TFILM,NDATA1,N0RDER,HFG) 

CALL  SEARCH(T,RH0,TI,NDATA2,N0RDER,RH0L) 

CALL  SEARCF(T,CP,TI,NCATA2,NCRDER,CL ) 

CALL  SEARCh(T,PRANTL,TI,NCATA2,NORDER,PRANTF ) 

CALL  SEARCH (TEMP, ALAMC A, TF ILM , NDATA1 , NORDER , ALAMF ) 

CCNVERT  ALAMF  TO  CALS/G-MOLE 
ALAMF=ALAMF*10. 

CCNVERT  RHCL  TC  GMS/CC 

RHCLB=RHCL*0.016 

CALCULATE  SURFACE  TENSION  FROM  WALDEN'S  EQUATION  (PERRY) 
IN  DYNES/CM 

ST=( ALAMF»RH0LB)/364. 

CCNVERT  SURFACE  TENSION  TO  LBS/FT 

ST=ST*68.5218E-06 

A=CL»DELTAT/HFG 

B=ST/RHCL 

FCIN3=(VISCF»FFG*A*A«A)/(  2  .  1  <5  7E-C6*DELTAT»SQRT  (  B  )  »PRANTF  «»5 

.1) 
RETURN 
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DELTAT  IS  NEGATIVE 


10 


FCINe=0. 

RETURN 

END 
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C  FCIF 


PURPOSE 
CALCU 
CCNVE 

USAGE 

FUNCT 

PARAKETE 
TBULK 
W 
TW1 

FUNCTION 
COI*NCN/C 

1PRANTLI8 
DATA  NOR 
CALL  SEA 
CALL  SEA 
CALL  SEA 
CALL  SEA 
CALL  SEA 
CALL  SEA 
FL0W=60. 
RE=4.856 
IF(RE.LE 

1  ( VISCF/V 
IF(RE.LE 

LPRANTF** 
IF(RE.GT 

KVISCF/V 
RETURN 
END 


LATES  THE  INSIDE  FILM  COEFFICIENT  DUE  TO  FORCED 
CTION  USING  THE  C I TTUS-BOELTER  EQUATION 

ICN  FCIF(TBULK,W,TW1  ) 
RS 

-  BULK  TEMPERATURE  OF  LIQUID,  DEG  F,  INPUT 

-  FLOW  RATE  OF  LIQUID,.  LBS/HR,  INPUT 

-  INNER  WALL  TEMPERATURE,.  DEG  F,  INPUT 

FCIF(TBULK, W, TW1 ) 
ATA2/NCATA2,T(8),RHO(8),CP(8) ,VISC(8) ,THCOND(8), 
)  ,EETA 
CER/4/ 

RCMT.RHO, TBULK, NDAT A2,NORDER , RHOF ) 
RCMT,VISC,TBULK,NDATA2,N0RDER,VISCF ) 
RCH(T»CP, TBULK,NCATA2,N0RDER,CPF ) 
RCH ( T, PR ANTLf TBULK, NDATA2fNORDERt PRANTF) 
RCF(T, THCCND,TBULK,NCATA2,NORDER,THCF) 
RCHT,VISC,TW1,NCATA2,.N0RDER,VISCW) 
*W 

7»FLOW/VISCF 

.2COO.  )FCIF=2.C9  3  414«THCF*RE«*0.545*PRANTF**0.4* 
ISCW)«*0.14 

.  5000..ANC.RE.GT.200C.  )  FC IF  =  0.8  8655* THCF*RE**0.66* 
:-.4*(  VISCF/V  ISCW)**0.  14 

.5000.  )FCIF  =  C.274  5000*THCF»RE**0.8*PRANTF**0.4* 
ISCW)  »»0. 14 
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C      PURPCSE 
C         CALCU 


C  FCO 

*«##**«»»»»«****«»#«*«#**«»*♦**##*»**#»»*»*»*»##*»# 

LATES  THE  OUTSIDE  (VAPOR  SIDE)  FILM  COEFFICIENT  USI 

NG 
USSELT  EQUATION 

ION  FC0(T0,TW2, IEFECT) 
RS 

-  MEAN  VAPCR  TEMPERATURE,.  DEG  F,  INPUT 

-  OUTER  WALL  TEMPERATURE,  DEG  F,  INPUT 
T  -  PARAMETER  INDICATING  THE  EFFECT 
AMS  REQUIRED 
UTINE  SEARCH 


10 


THE  N 

USAGE 

FUNCT 

PARAMET6 
TO 
TW2 
IEFEC 

SUBPRCGR 
SUERC 

FUNCTION 
COVMCK/C 

1PRANTLU 
CONPCN/C 

1ALAMDA(2 
DATA  NOR 
CALL  SEA 
TFIL^=TC 
CALL  SEA 
CALL  SEA 
CALL  SEA 
DELTAT=T 
AL=9.5 
IFUEFEC 
G=4.17E0 
IF(DELTA 
FC0=0.1E 
RETURN 
FC0=0.9A 
RETURN 
END 


FCO(T 
ATA2/N 
),BETA 
ATA1/N 
2), PRE 
CER/4/ 
RCF(TE 
-0.75* 
RCH(T, 
RCHtT, 
RCH(T, 
C-TW2 


0,TW2,  IEFECT) 

DATA2,T(8),RH0(  8  )  ,  CP  (  8  )  ,.VI  SC  (  8  )  ,THC0ND(8)  , 

CAT  A 1, TEMP (22) , THALP ( 22 ) , THALP V ( 22 ) , 
SS(22),RHOV(22) 

MP,  ALAMCA,T0,NCATA1,.N0RDER,ALAMF) 
(T0-TW2) 

THC0ND,TFILM,NCATA2,N0RDER,THCF) 
RH0,TFILM,NDATA2,N0RDER,RH0F) 
VISC,TFILM,NDATA2,N0RDER,VISCF ) 


T.EQ.2)AL=2.3125 

8 

T.GT.O.IGO  TC  10 

+  40 

3*(THCF**3*RF0F*»2*G*ALAMF/(DELTAT*AL*VISCF))**0.25 
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U 


c 

c  *  »  »*  »  * 
C      P 

c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 


FCIN 

URPCSE 

CALCULATES  THE  INSIDE  FILM  COEFFICIENT  FOR  NATURAL 
CCNVECTICN  USING  THE  NATURAL  CONVECTION  EQUATION 

SAGE 

FUNCTICN  FCIN(TBULK,TW,H,  IEFECT) 

ARAMETERS 

T8ULK   -  BULK  TEMPERATURE  OF  LIQUID,  DEG  F,  INPUT 
TW      -  INNER  WALL  TEMPERATURE,  DEG  F,  INPUT 
H       -  HOLD-UP,  LBS,  INPUT 
IEFECT  -  PARAMETER  INDICATING  THE  EFFECT 

UBPRCGRANS  REQUIRED 
SUBROUTINE  SEARCH 
FUNCTICN  HEIGHT 


F 
C 
1A 
D 
T 
I 
C 

c 
c 
c 

X 

I 

G 

Y 
I 
1* 
I 
I 
R 

F 
R 
F 
R 
F 
E 
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UNCTICN 

CNMCN/L 

LAMCA(2 

ATA  NOR 

FILM=0. 

F(  (Tk-T 

ALL  SEA 

ALL  SEA 

ALL  SEA 

ALL  SEA 

L=HEIGh 

F(A3S(X 

R=(XL»* 

=GR*PRA 

F( (Y.GE 

(Y»*0.2 

FtY.GE. 

FIY.LT. 

ETURN 

RITE(6, 

CIN=(IC 

ETURN 

CIN=  . 

ETURN 

0RMAT(5 

ND 


FCIN( 
ATA1/N 
2)  ,PRE 
DERM/ 
5MTBU 
BULK) . 
RCHIT, 
RCH(T, 
RCH(T, 
RCH(T, 
T(1,H, 
D.LT. 
3*RH0F 
NTF 

.1.E4) 
5) 

3.5E7) 
1.E4JG 


TBULK,TW,H, IEFECT) 

CAT Al, T EMP { 22 ),THALP( 22 ),THALPV< 22) , 

SSI 22),RH0V(22  ) 

LK+TW) 

LE.O.)GO  TO  10 

THC0NC,TFILM,NCATA2,NORDER,THCF) 

RH0,TFILM,NDATA2,N0RDER,RH0F) 

VISC,TFILM,NCATA2,N0RDER,VISCF) 

PRANTL,TFILM,NCATA2,N0RDER,PRANTF) 

TBULK,TW,0.,0.,  IEFECT) 

•1E-70JG0  TO  10 

*RH0F*4.17E8*BETA*(TW-TBULK) ) / ( VI SCF»VI SCF ) 

.AMU.(Y.LT.3.5E7))FCIN=( (0.5  5»THCF )/XL) 

FCIN=(  (0.13*THCF)/XL)»(Y»*C33  3) 
0  TO  2C00 


2100) 

.55*THCF)/XL  )*(Y»»C25) 


X.'Y  IS  OUT  CF  RANGE'  ) 
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C*«««**«****«* 


£«#*«#«*»»**»* 


10 


20 


THIS 

DESIR 
USAGE 

CALL 
PARAI^ETE 

X 

N 

Y 

IPCSI 


M 
Z 

SUBRCUTI 

CINENSIC 

DO  10  1= 

Z(I)=X( I 

CONTINUE 

IC0UNT=1 

J=IPCSIT 

Z( J)=Y( I 

ICOUNT=I 

IFdCCUN 

RETURN 

END 


MAPI 

SUBROUTINE  MAPS  THE  ELEMENTS  OF  A  VECTOR  INTO 
EC  POSITIONS  IN  ANOTHER  VECTOR 

MAPK  X,N, Yf IPOSIT,M,Z) 
R  S 

-  INTO  VECTOR  INTO  WHICH  ELEMENTS  ARE  MAPPED 

-  DIMENSION  OF  X  AND  Z,  INPUT 

-  INPUT  VECTOR  TO  REPLACE  ELEMENTS  IN  X 

T  -  VECTOR  CONTAINING  INTEGERS  INDICATING  WHICH 
ELEMENTS  IN  X  ARE  TO  BE  REPLACED  BY  THE 
ELEMENTS  OF  Y 

-  DIMENSION  OF  Y  AND  IPOSIT,  INPUT 

-  OUTPUT  VECTOR 


NE  MAPI  (X,N,Y,  IPOSIT, M, Z  ) 
N  X(1),Y(  1),  IPOSITU  ),Z(  1) 
1,N 
) 


(  ICOUNT) 
COUNT  ) 
COUNT  41 
T.LE.M)GO  TO  20 
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